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SUMMARY 

This  report  describes  the  results  of  the  second  year  of  a 
three-year  research  program  to  Investigate  methods  for  obtaining 
diffraction-limited  images  of  space  objects,  despite  the  turbulent 
atmosphere,  by  reconstructing  Images  from  data  provided  by  optical 
Interferometers  (particularly  stellar  speckle  interferometry). 
Accomplishments  Include  the  following.  (1)  Improved  image 
reconstruction  algorithms  were  developed.  (2)  A  better  understanding  of 
modes  of  stagnation  of  algorithms  was  developed.  (3)  The  performances 
of  the  shlft-and-add  Image  formation  method  and  of  one  recursive 
algorithm  were  Investigated.  (4)  A  second  recursive  algorithm  was  shown 
to  suffer  from  a  uniqueness  problem.  (5)  A  potential  new  remote  sensing 
application  of  the  Iterative  reconstruction  algorithms  was  explored. 
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DIFFRACTION-LIMITED  IMAGING  OF  SPACE  OBJECTS  II 
1  March  1983  to  29  February  1984 


INTRODUCTION  AND  OBJECTIVES 

This  report  describes  the  results  of  the  second  year’s  effort  in  a 
three-year  research  program  to  Investigate  methods  of  obtaining 
diffraction-limited  Images  of  space  objects,  despite  the  turbulent 
atmosphere,  by  reconstructing  Images  from  data  provided  by  optical 
Interferometers  (particularly  stellar  speckle  interferometry). 

Atmospheric  turbulence  typically  limits  the  angular  resolution  of 
earth-bound  optical  telescopes  to  one  second  of  arc  or  worse,  which  is 
fifty  times  poorer  than  the  theoretical  diffraction  limit  of  a  5-meter 
optical  telescope.  It  is  possible  to  gather  diffraction-limited 
information  through  the  turbulent  atmosphere  by  a  variety  of 
interferometric  techniques,  including  Michelson  stellar  interferometry 
[1],  intensity  interferometry  [2],  amplitude  interferometry  [3],  and 
stellar  speckle  interferometry  [4,  5].  However,  this  diffraction- 
limited  information  is  in  the  form  of  the  modulus  (magnitude)  of  the 
Fourier  transform  of  the  object  being  viewed.  Until  recently  only  the 
autocorrelation  of  the  object,  but  not  the  object  itself,  could  be 
reconstructed  from  this  data,  except  for  special  cases. 


In  recent  years  an  iterative  method  [6-9]  has  been  developed  for 
reconstructing  an  object  from  its  Fourier  modulus,  thereby  making 
possible  the  reconstruction  of  diffraction-limited  imagery  from 
interferometer  data.  The  algorithm  utilizes  the  measured  Fourier 
modulus  data  as  well  as  (1)  the  a  priori  information  that  the  object's 
spatial  (or  angular)  brightness  distribution  is  a  nonnegative  function 


and  (2)  Information  about  the  object's  diameter  which  can  be  computed 
from  the  autocorrelation  function.  The  algorithm  and  Its  numerous 
applications  Is  described  In  detail  In  Appendix  A  [9]. 

The  goal  of  the  program  is  to  further  investigate  and  develop  this 
method  of  obtaining  diffraction-limited  Images.  Included  in  the 
three-year  program  are  Investigations  Into  improving  the  reconstruction 
algorithm,  developing  methods  for  processing  noisy  astronomical  data, 
studying  the  uniqueness  of  the  reconstruction,  and  investigating  ways  to 
Increase  the  spectral  bandwidth  of  stellar  speckle  interferometry.  In 
the  second  year  of  the  effort,  the  emphasis  was  on  developing  new  and 
Improved  reconstruction  algorithms.  Initial  studies  of  the  uniqueness 
problem  and  of  the  properties  of  astronomical  data  were  also  begun. 

The  research  accomplishments  for  the  second  year  are  summarized  in 
Section  2  and  are  described  In  more  detail  In  Sections  3  through  8  and 
In  the  Appendices.  References  are  listed  In  Section  9. 
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RESEARCH  ACCOMPLISHMENTS 

The  second  year  of  research  effort  can  be  divided  Into  seven  major 
topics. 


1.  The  new  recursive  algorithm  [10]  described  in  last 
year's  report  [11]  was  implemented  and  tested  both  on 
noise-free  and  on  noisy  Fourier  modulus  data. 

2.  Improvements  in  the  iterative  Fourier  transform 
reconstruction  algorithm  [6-8]  were  made  that  enable  one  to 
reconstruct  difficult  objects  that  previously  defied 
reconstruction  attempts. 

3.  Alternative  Iterative  algorithms  were  devised. 

A.  Investigations  were  made  into  the  problem  of  stripes 
In  the  reconstructed  Images. 

5.  The  shift-and-add  algorithm  was  implemented  and 
tested  on  a  complicated  extended  object. 

6.  Results  were  obtained  indicating  a  possible  new 
remote  sensing  application  of  the  iterative  reconstruction 
algorithm. 

7.  Recently  published  claims  regarding  the  uniqueness  of 
phase  retrieval  when  the  edges  of  the  object  are  known  were 
shown  to  be  false  by  counterexample  [12], 

Recent  publications  arising  from  this  work  are  References  10  and 
12-18.  Reference  9  is  noted  as  a  recent  related  publication  arising 
from  a  previous  research  program  [19]  and  is  included  as  Appendix  A. 
References  10,  16,  17,  18  and  12  are  included  as  Appendices  B,  C,  D,  E 
and  F,  respectively. 
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The  seven  topics  are  briefly  described  in  the  remainder  of  this 
section  and  are  described  in  more  detail  In  Sections  3  to  8  and  in  the 
Appendices. 

2.1  NEW  RECURSIVE  ALGORITHM 

As  described  in  last  year's  report  [11]  and  in  Appendix  B,  a  new 
recursive  algorithm  has  been  developed  which  is  capable  of 
reconstructing  an  object  from  Its  autocorrelation  function,  which  can  be 
computed  from  the  modulus  of  Its  Fourier  transform.  It  works  for 
objects  having  latent  reference  po1nts--unresol ved  points  within  the 
object  field  that  are  not  sufficiently  far  from  the  main  part  of  the 
object  to  satisfy  the  condition  for  holography,  but  sat4  ;  weaker 
conditions.  The  recursive  algorithm  was  coded  on  ymputer  and 
exercised  on  two  different  types  of  objects  using  autocorrelations 
having  a  variety  of  signal -to-nol se  ratios.  As  expected  „..c  recursive 
algorithm  was  fairly  sensitive  to  noise,  making  It  less  practical  for 
real-world  applications  than  the  Iterative  Fourier  transform  algorithm 
[6-8].  Improvements  were  made  In  the  recursive  algorithm  to  make  it 
somewhat  less  sensitive  to  noise.  A  more  detailed  description  of  this 
work  is  given  In  Section  3. 

2.2  IMPROVEMENTS  IN  THE  ITERATIVE  ALGORITHM 

The  Iterative  Fourier  transform  algorithm  has  been  particularly 
successful  on  objects  having  complicated  shapes,  such  as  satellites  [6, 
7,  17].  However,  for  some  types  of  objects,  such  as  those  whose  support 
fills  a  square,  the  algorithm  has  a  tendency  to  stagnate  without  finding 
a  solution  [20].  Two  modifications  of  the  basic  approach  were 
demonstrated  for  overcoming  this  problem.  The  first  was  to  employ  the 
defogglng  method  of  Bates  and  Fright  [21].  The  second  modification  Is 
to  break  the  symmetry  of  the  partially  reconstructed  image  in  order  to 
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allow  the  algorithm  to  converge  to  one  of  the  two  possible  solutions. 
This  is  described  further  and  an  example  Is  shown  In  Section  4. 

2.3  ALTERNATIVE  ITERATIVE  ALGORITHM 

The  theoretical  justification  for  the  input-output  iterative 
Fourier  transform  algorithm  [6-9]  alludes  to  a  control  theory 
poi nt-of-view.  Yet  rigorous  control  theory  had  not  actually  been 
applied  to  the  problem.  Alternative  algorithms  based  on  control  theory 
are  presented  in  Section  5.  Further  elaboration  of  these  algorithms  and 
their  implementation  and  testing  will  be  required  to  determine  whether 
they  will  offer  improved  performance  over  existing  algorithms. 

2.4  INVESTIGATION  OF  STRIPES 

In  some  cases  the  output  of  the  Iterative  algorithm  has  the 
appearance  of  original  object  but  with  a  pattern  of  low-contrast  stripes 
superimposed  [22,  17].  The  phenomenon  of  stripes  appearing  in  the 
reconstructed  image  was  extensively  investigated.  Properties  of  the 
phase  of  the  Fourier  transform  of  the  striped  image  were  studied. 
Several  methods  for  removing  the  stripes  were  investigated.  We  feel 
that  we  are  on  the  threshold  of  solving  this  problem,  as  described  in 
Section  6.  When  this  problem  is  completely  solved,  then  it  will  be 
possible  to  answer  the  question  of  the  uniqueness  of  the  reconstructed 
image  (see  Appendix  D)  more  definitively. 

2.5  SHIFT-AND-ADD 

The  shi f t-and-add  method  of  imaging  from  short-exposure 
astronomical  Images  has  In  the  past  been  exercised  primarily  for  very 
simple  objects  having  In  their  f leld-of-view  very  bright  unresolved 
points  [23,  24].  The  shif t-and-add  method  was  attempted  both  on  a 


Terim 


point-like  object  and  on  a  more  realistic  object— a  satellite  having 
strong  glints.  Although  the  result  for  the  point-like  object  was  very 
good,  the  result  for  the  extended  object  was  very  poor,  indicating  that 
the  shift-and-add  method  Is  not  appropriate  for  complicated  extended 
objects.  These  results  are  shown  in  Section  7. 

2.6  NEW  ITERATIVE  APPLICATION 

A  new  remote  sensing  application  for  the  iterative  Fourier 
transform  algorithm  was  developed  under  ERIM  internal  funding  [25].  It 
permits  the  operation  of,  say,  a  synthetic  aperture  radar  system  having 
reduced  performance  requirements  for  the  phase  stability  of  its  local 
oscillator  and  motion  compensation.  It  might  also  be  useful  for  the 
electron  microscopy  phase  retrieval  problem.  It  involves  the  iterative 
retrieval  of  phase  using  a  single  intensity  measurement  plus  a  shape 
constraint  on  the  object  or  upon  the  pattern  of  radiation  by  which  the 
object  is  Illuminated.  Under  the  present  program  the  issue  of  the  shape 
constraint  was  explored  further.  It  was  found  that  certain  Interesting 
shapes  are  sufficient  for  reconstructing  a  complex-valued  object 
function  from  the  magnitude  of  its  Fourier  transform .  The 
reconstruction  algorithm  and  some  reconstruction  results  are  shown  in 
Section  8. 

2.7  AMBIGUITY  OF  PHASE  RETRIEVAL  USING  BOUNDARY  CONDITIONS 

Claims  have  been  made  that  an  object  can  be  uniquely  reconstructed 
from  its  Fourier  modulus  via  the  autocorrelation  function  if  the  values 
of  the  edges  of  the  object  are  known  [26].  It  is  shown  In  Appendix  F 
that  knowledge  of  the  autocorrelation  function  and  of  the  boundary 
values  of  the  object  are  not  sufficient  to  uniquely  specify  the  object 
In  all  cases.  It  Is  further  shown  how  and  why  the  recursive  algorithm 
of  Hayes  and  Quatieri  [26]  fails  for  the  nonunique  cases.  Recently  it 
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EXPERIMENTAL  RESULTS  USING  NEW  RECURSIVE  ALGORITHM 

As  described  In  Appendix  B,  a  new  recursive  algorithm  has  been 
developed  for  reconstructing  an  object  from  Its  autocorrelation 
function,  which  can  be  computed  from  the  modulus  of  the  Fourier 
transform  of  the  object.  It  Is  applicable  to  objects  having  latent 
reference  points  [10],  and  knowledge  of  the  support  of  the  object  may  be 
required.  In  this  section  examples  of  reconstruction  experiments  using 
the  recursive  algorithm  are  shown. 

Figure  3-1  shows  results  of  the  recursive  reconstruction  algorithm 
for  which  the  Fourier  modulus  (or  autocorrelation)  data  was  corrupted 
with  varying  amounts  of  noise.  The  object  consists  of  an  equilateral 
right  triangle  of  16  pixels  on  each  side  having  a  brighter  rectangle  and 
a  brighter  square  Imbedded  In  It.  It  Is  assumed  known  that  the  object's 
support  Is  the  triangle.  Figure  3-l(a)  shows  the  original  object. 
Figure  3- 1 ( b ) ,  (c)  and  (d)  are  the  reconstructed  Images  when  the 
root-mean-squared  (RMS)  error  of  the  Fourier  modulus  data  was  0.005175, 
0.05585  and  0.01795,  respectively.  The  RMS  error  of  these  reconstructed 
Images  is  0.0400,  0.6088  and  0.1390,  respectively.  That  Is,  for  the 
noisiest  case  of  Figure  3-l(c),  a  5.6  percent  error  in  the  data  resulted 
In  a  60.9  percent  error  In  the  reconstructed  Image. 

To  get  a  feel  for  how  bad  a  60  percent  error  Is,  consider  the 
following.  Suppose  the  object  were  constant,  equal  to  unity,  over  the 
known  region  of  support  (within  the  triangle).  If  the  reconstructed 
Image  were  a  set  of  random  numbers  uniformly  distributed  between  0  and 
b,  then  the  rms  error  for  the  optimum  value  of  b  can  be  shown  to  be  50 
percent.  That  Is,  the  reconstructed  Image  shown  In  Figure  3- 1 ( c ) , 
having  RMS  error  of  60  percent.  Is  worse  than  a  reconstructed  Image 
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Figure  3-1.  Images  Reconstructed  Using  Recursive  Algorithm. 

(a)  Original  object;  (b)-(d)  images  reconstructed  from 
autocorrelations  having  root-mean-squared  error 

(b)  0.0400,  (c)  0.6088,  and  (d)  0.01795. 


consisting  of  random  numbers. 

By  comparison.  In  experiments  using  the  Iterative  Fourier  transform 
algorithm  (on  a  different  object),  a  5*  error  In  the  data  resulted  In  a 
20*  error  In  the  reconstructed  Image  [22].  Therefore  from  this  limited 
experience  It  appears  that  the  recursive  algorithms  Is,  as  predicted 
[10],  highly  sensitive  to  noise.  The  Iterative  Fourier  transform 
algorithm  would  appear  to  be  the  preferred  method  of  Image 
reconstruction. 

For  the  triangular  support  case,  one  can  generate  as  many  as  three 
separate  estimates  for  each  value,  one  associated  with  each  of  the  three 
corner  pixels.  In  the  absence  of  noise  these  three  reconstructions  are 
Identical,  but  with  noise  present  they  will  In  general  be  different,  and 
an  Improved  algorithm  would  decrease  noise  effects  by  averaging  the 
three  estimates.  This  method  was  tried  and  the  following  was 
discovered.  The  computation  of  a  value  depends  not  only  the  corner 
pixel  but  also  on  a  number  of  previously  reconstructed  values  as  well. 
In  the  presence  of  noise  each  of  these  reconstructed  values  will  have 
some  error  associated  with  It,  and  the  more  of  them  that  are  used  to 
reconstruct  a  new  value,  the  more  error  that  new  value  will  have.  It 
was  found  that  the  estimate  that  Is  generated  from  the  maximum  number  of 
previously  reconstructed  values  has  accumulated  so  much  error  that 
Including  Its  value  In  the  average  degrades  rather  than  Improves  the 
reconstruction.  The  optimum  number  of  estimates  to  use  was  found  to 
depend  on  the  slgnal-to-nolse  ratio  and  on  distance  from  the  edges  of 
the  triangle.  For  most  values,  only  one  or  two  estimates  was  optimal. 
This  resulted  In  a  modest  Improvement  over  the  algorithm  employing  only 
a  single  estimate. 


4 

IMPROVEMENTS  IN  THE  ITERATIVE  FOURIER  TRANSFORM  ALGORITHM 


Although  the  iterative  Fourier  transform  algorithm  has  been  shown 
to  be  successful  for  space  objects  such  as  satellites  [6,  7,  17],  It  can 
have  problems  converging  for  some  other  types  of  objects.  For  example, 
for  the  object  shown  In  Figure  4-l(a),  a  picture  of  a  bird  bounded  by  a 
square,  the  algorithm  has  a  strong  tendency  to  stagnate  without  finding 
a  solution  [20].  Two  methods  were  demonstrated  for  overcoming  this 
problem  and  converging  to  a  solution:  the  defogglng  method  of  Bates  and 
Fright  [21]  and  a  new  method  of  temporarily  using  a  reduced-area 
asymmetric  support  constraint.  The  latter  method  seems  to  be  the  more 
Important  of  the  two. 

The  defogglng  method  attempts  to  compensate  for  the  fact  that  a 
low-contrast  object  on  a  bright  background  causes  very  little 
"Interference,"  that  Is,  its  Fourier  transform  has  not  much  structure. 
The  defogglng  method  consists  of  reducing  the  large  central  lobe  of  the 
Fourier  modules,  raising  the  relative  values  at  the  higher  spatial 
frequencies.  In  the  image  domain  this  corresponds  to  reducing  any 
slowly-varying  bias-like  (or  fog)  component  of  the  Image,  thereby 
emphasizing  the  f Iner-structure  details.  Phase  retrieval  algorithms 
that  work  poorly  on  a  low-contrast  image  tend  to  work  better  on  the 
defogged  version  of  the  Image.  After  the  defogged  Image  Is 
reconstructed,  the  slowly- varying  fog  component  Is  added  back  In.  As  a 
final  step  the  refogged  image  Is  refined  by  further  iterations  of  phase 
retrieval . 

A  reduced-area  asymmetric  support  constraint  Is  used  for  types  of 
objects  which  often  cause  the  Iterative  Fourier  transform  algorithm  to 
stagnate  on  a  partially  reconstructed  Image.  Recall  that  there  Is 


Figure  4-1.  Use  of  Asymmetric  Support  Constraint,  (a)  Original 

object;  (b)  output  image  from  iterative  Fourier  transform 
algorithm  which  has  stagnated;  (c)  output  image  after 
application  of  reduced-area  support  constraint  followed  by 


always  a  two-fold  ambiguity:  f(-x,  -y)  and  f(x,  y)  have  the  same  Fourier 
modulus.  The  ambiguous  Image  f(-x,  -y)  Is  just  f(x,  y)  rotated  by  180° 
which  Is  equivalent  to  being  reflected  through  the  origin.  When  there 
Is  a  symmetric  support  as  In  the  case  of  the  object  shown  In  Figure 
4-l(a),  the  algorithm  may  stagnate  with  an  output  such  as  the  one  shown 
In  Figure  4-l(b),  which  has  features  of  both  the  object  and  the  180° 
rotated  object.  The  output  Image  changes  little  with  further 
Iterations.  Apparently  the  algorithm  gets  stuck  when  It  Is  half-way 
between  the  two  different  solutions:  It  is  unable  to  shake  one  off  and 
converge  to  the  other.  In  the  case  where  an  asymmetric  support 
constraint  Is  known,  this  particular  mode  of  stagnation  tends  not  to 
occur  since  the  asymmetric  support  constraints  moves  the  solution  toward 
f(x,  y)  and  away  from  f(-x,  -y). 

The  method  of  using  a  reduced-area  asymmetric  support  constraint  is 
as  follows.  A  support  constraint  Is  defined  that  Includes  one  side 
(Including  edges)  of  the  object  but  not  the  other  side  and  edges.  This 
support  constraint  Is  chosen  to  be  smaller  than  the  known  support  of  the 
object  and  to  be  as  asymmetric  as  possible,  so  that  It  has  little  In 
common  with  the  180°-rotated  version  of  the  support  constraint.  A  few 
Iterations  are  then  performed  with  the  reduced-area  asymmetric  support 
constraint  (rather  than  using  the  correct  support  constraint).  It  Is 
hoped  that  this  causes  one  of  the  two  Images,  f(x,  y)  or  f(-x,  -y)  to  be 
preferentially  enhanced  over  the  other.  After  switching  back  to  the 
correct  support  constraint,  either  f(x,  y)  or  f(-x,  -y)  will  be  strong 
enough  compared  with  the  other  that  upon  further  iterations  the 
algorithm  converges  to  the  stronger  one  and  away  from  the  weaker  one. 

Figure  4- 1 ( c )  shows  the  reconstructed  output  Image  after  using  both 
the  defogglng  method  and  the  reduced-area  asymmetric  support  constraint 
for  a  few  Iterations  then  continuing  with  further  Iterations  using  the 
correct  support  constraint  and  the  original  Fourier  modulus  data. 


Comparing  It  with  the  output  Image  shown  In  Figure  4-l(b),  It  Is  seen 
that  these  techniques  yielded  much  better  results  In  this  case. 

These  methods  have  been  exercised  In  only  very  limited  circum¬ 
stances  and  have  not  yet  been  optimized  and  automated.  Further  work  to 
develop  these  promising  methods  Is  clearly  called  for.  Such  auxiliary 
procedures  are  not  necessary  for  the  objects  that  are  easier  to 
reconstruct  but  are  very  Important  for  the  more  difficult  cases. 


5 

ALTERNATIVE  ITERATIVE  ALGORITHMS 


The  iterative  Fourier  transform  algorithm,  which  is  described  in 
detail  In  Appendix  A,  works  very  well  in  a  wide  range  of  situations  but 
converges  slowly  or  not  at  all  In  some  cases.  Furthermore  it  is  always 
desired  to  arrive  at  a  solution  using  fewer  Iterations  and  less  computer 
time.  For  these  reasons  we  are  always  looking  for  ways  to  Improve  the 
existing  algorithms  or  devise  alternative  algorithms  that  converge 
faster.  The  two  algorithms  shown  In  Figures  5-1  and  5-2  are  examples  of 
alternative  algorithms  that  have  been  conceived.  They  were  arrived  at 
from  the  point  of  view  of  control  theory. 

In  the  first  algorithm,  depicted  In  Figure  5-1,  it  Is  assumed  that 
the  individual  sldelobes  of  the  complex  Fourier  transform  of  the  object 
can  be  modelled  by  a  fairly  simple  mathematical  formula  having  a  small 
number  of  free  parameters .  By  curve  fitting  each  lobe  of  the  Fourier 
modulus  (amplitude)  to  the  model,  one  could  determine  the  parameters  and 
thereby  determine  the  phase.  One  would  first  curve  fit  one  larger  lobe, 
compute  the  magnitude  of  the  model  from  the  fitted  parameters,  and 
subtract  that  model  from  the  modulus  measurement.  Smaller  lobes  would 
be  curve-fitted  and  subtracted  resurslvely  from  the  residual  modulus. 
After  all  the  lobes  are  modelled,  the  corresponding  phase  would  be 
combined  with  the  measured  modulus  and  the  Image  would  be  computed  by 
Inverse  Fourier  transformation.  It  Is  yet  to  be  determined  whether  the 
Fourier  transform  can  be  modelled  as  described  above. 

In  the  second  algorithm,  shown  In  Figure  5-2,  the  difference  in  the 
phase  of  the  Fourier  transform  of  the  current  estimate  and  that  of  the 
previous  estimate  Is  multiplied  by  a  gain  factor  (K)  and  added  to  the 
previous  phase  estimate.  This  Is  similar  to  previous  iterative 
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Figure  5-2.  Alternative  Reconstruction  Algorithm  -  ] 


algorithms  except  that  the  roles  of  the  two  domains  are  reversed. 

Both  the  methods  described  above,  as  well  as  others,  merit  further 
research  and  implementation. 
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INVESTIGATION  OF  THE  STRIPES  PHENOMENON 

In  a  number  of  cases  the  Iterative  Fourier  transform  algorithm  has 
converged  almost  all  the  way  to  a  solution,  but  then  stagnates  at  an 
output  Image  that  looks  Tike  the  original  object  but  having  a  set  of 
stripes  superimposed  [17,  22]  (see  Appendix  D).  In  most  cases  the 
stripes  are  of  such  low  contrast  as  to  be  hardly  noticeable,  but 
occasionally  the  contrast  of  the  stripes  has  been  high  enough  to  be 
objectionable.  Since  the  Fourier  transform  pair  is  not  in  perfect 
agreement  with  the  data  and  constraints  in  this  condition,  the  striped 
images  Is  at  a  local,  rather  than  the  global,  minimum  of  the  error,  and 
therefore  it  does  not  represent  a  lack  of  uniqueness.  Although  earlier 
attempts  at  solving  this  problem  (the  stagnation  at  a  striped  imaged) 
failed,  we  are  currently  developing  methods  that  will  eliminate  the 
stripes. 

Since  an  output  Image  having  stripes  has  a  Fourier  modulus  equal  to 
the  measured  (assumed  to  be  the  correct)  Fourier  modulus,  the  stripes 
must  be  due  to  the  effects  of  phase  errors.  The  phase  errors  must  be 
located  in  small  regions  of  the  Fourier  domain  in  order  to  produce  such 
a  regular  striped  pattern,  with  the  locations  of  the  regions  in  the 
Fourier  domain  being  related  to  the  spatial  frequency  (spacing)  and 
orientation  (angle)  of  the  stripes. 

A  first  attempt  at  eliminating  the  stripes  was  to  add  noise  to  the 
Input  Image  after  stagnation  at  a  striped  output  had  occurred.  The  hope 
was  that  the  added  noise  would  move  the  solution  far  enough  away  from 
the  local  minimum  so  that  further  Iterations  of  the  Iterative  Fourier 
transform  algorithm  would  lead  to  the  global  minimum  rather  than  falling 
back  into  the  same  local  minimum.  When  this  was  tried  it  was  found  that 
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after  further  Iterations  the  algorithm  did  Indeed  fall  back  Into  the 
same  local  minimum,  even  when  the  amount  of  noise  added  was  very  large. 


A  second  attempt  at  solving  the  stripes  problem  relied  on  the 
knowledge  that  the  phase  error  was  a  localized  phase  error  In  the 
Fourier  domain.  It  was  found  that  If  a  constant  phase  was  added  to  the 
phase  of  the  Fourier  transform  of  the  object  in  a  given  region  of  the 
Fourier  domain  (and  In  order  to  preserve  the  Hermitian  property  of  the 
Fourier  transform,  the  same  constant  phase  was  subtracted  in  the 
symmetric  region  of  the  Fourier  domain),  then  the  corresponding  image 
would  look  like  the  original  object  but  with  a  pattern  of  stripes 
superimposed.  The  resulting  synthesized  Images  thus  produced  had  an 
appearance  very  much  the  same  as  the  striped  images  produced  by  the 
Iterative  algorithm.  However,  when  these  Images  were  used  as  the  Input 
to  the  iterative  algorithm,  the  synthesized  stripes  immediately  went 
away  and  the  algorithm  quickly  converged  to  the  true  image.  This 
constrasted  sharply  with  the  stripes  produced  by  the  iterative 
algorithm,  which  would  not  go  away.  Therefore  the  phase  errors  that 
cause  the  stripes  problem  are  more  complicated  than  simple  constant 
phase  errors  over  some  region  of  the  Fourier  domain.  This  was  further 
shown  by  attempts  to  eliminate  the  stripes  by  adding  various  constant 
phase  errors  at  appropriate  regions  of  the  Fourier  domain.  All  such 
attempts  failed  to  remove  the  stripes. 


A  third  attempt  at  removing  the  stripes  Involved  the  addition  of 
Blaschke-like  phase  functions  to  the  Fourier  transform.  A  Blaschke-like 
phase  Is  the  phase  of  the  unity-modulus  function 
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where  u  Is  a  Fourler-domaln  coordinate  and  z  =  a  +  1b.  The  Blaschke- 
1 1 ke  phase  Is  global  In  effect  but  has  Its  most  rapid  variation  within  a 
region  about  u  «  b.  Varying  phase  error  corrections  of  this  form  were 
also  found  to  be  unsuccessful  In  solving  the  stripes  problem. 

Most  recently  we  have  developed  two  procedures  that  should  solve 
the  stripes  problem  In  most  cases.  They  are  based  upon  two  facts: 

1.  The  phase  errors  that  produce  the  stripes  are  located 
In  small  regions  of  the  Fourier  domain,  and 

2.  Output  Images  arrived  at  by  the  Iterative  Fourier 
transform  algorithm  started  with  different  Initial  Inputs  of 
random  numbers  are  unlikely  to  have  the  same  pattern  of 
stripes. 

In  the  first  method,  three  output  Images  are  produced  by  the 
Iterative  Fourier  transform  algorithm  using  three  different  initial 
inputs  of  random  numbers.  The  three  output  images  are  translated  so  as 
to  be  centered  at  the  same  point  In  order  to  remove  any  linear  phase 
difference  In  their  Fourier  transforms.  Then  at  each  point  in  the 
Fourier  domain,  the  complex  values  of  the  three  Fourier  transforms  are 
compared.  The  value  whose  distance  from  the  other  two  values  Is  the 
greatest  Is  discarded  and  a  new  value  Is  formed  by  taking  the  average  of 
the  remaining  two  (closest)  values.  In  this  manner.  If  In  a  given 
region  of  the  Fourier  domain  one  of  the  three  has  a  phase  error  (related 
to  the  stripes  or  otherwise),  that  phase  error  Is  eliminated. 

In  the  second  method  only  two  different  output  images  need  to  be 
produced  by  the  Iterative  Fourier  transform  algorithm.  Although  the 
stripes  are  typically  of  highest  contrast  where  the  object  Is  brightest, 
they  also  exist  outside  the  known  support  of  the  object.  So  by  Fourier 
transforming  the  region  of  the  output  image  having  only  stripes  (outside 
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the  support  of  the  object),  the  regions  of  the  Fourier  domain  having  the 
phase  error  can  be  Identified.  Then  a  new  phase  estimate  Is  made  by 
using  the  phase  of  the  Fourier  transform  of  the  first  output  Images 
where  it  Is  not  Influenced  by  the  phase  error,  and  using  the  phase  of 
the  Fourier  transform  of  the  second  output  Image  where  the  first  was 
Influenced  by  the  phase  error. 

In  both  methods  above,  after  the  new  estimate  Is  formed,  further 
iterations  should  be  performed  to  allow  It  to  converge  closer  to  a 
solution. 

These  methods  of  correcting  the  stripes  are  automatic  In  the  sense 
that  no  human  judgement  or  decisions  are  required  during  their 
operation. 

Both  methods  were  exercised  on  a  single  example  and  were  found  to 
perform  very  well.  Further  experimentation  with  these  methods  is 
required  to  determine  their  effectiveness  in  a  wider  variety  of 
circumstances. 
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EXPERIMENTAL  RESULTS  USING  SHIFT-AND-ADD 


Shlft-and-add  [23,  24]  Is  a  method  of  reconstructing  Images  of 
astronomical  objects  from  multiple  short-exposure  images.  It  consists 
of  shifting  all  the  Images  so  that  their  maximum  values  all  lie  at  the 
same  coordinate,  then  adding  (or  averaging)  them  all.  This  has  been 
shown  to  work  well  for  objects  consisting  of  a  collection  of  delta 
functions  (points)  [23,  24],  but  It  was  not  demonstrated  for  realistic 
extended  objects,  such  as  satellites.  We  Implemented  the  shlft-and-add 
method  on  the  computer  and  exercised  It  both  on  an  object  consisting  of 
a  collection  of  delta  functions  and  on  an  extended  object. 

Figure  7 - 1 ( a )  shows  an  object  consisting  of  three 
delta-function-like  points  having  relative  brightness  of  100:20:10. 
Simulated  point-spread  functions  from  a  telescope  Including  atmospheric 
turbulence  [22]  were  convolved  with  the  object  to  arrive  at  simulated 
blurred  Images,  an  example  of  which  Is  shown  In  Figure  7-l(b).  Figure 
7- 1(c)  shows  the  result  of  shifting  and  adding  156  blurred  Images.  The 
three  points  can  be  clearly  seen,  although  they  are  slightly  blurred  and 
each  Is  surrounded  by  a  fog.  We  took  shlft-and-add  one  step  further  by 
combining  it  with  a  form  of  subtractive  deconvolution  related  to  the 
CLEAN  method.  An  estimate  of  the  effective  point-spread  function  was 
found  by  applying  shlft-and-add  to  a  single  point.  The  outputs  from 
shift  and  add  were  then  CLEANed  by  subtracting  from  the  Image  a  version 
of  this  point-spread  function  shifted  and  scaled  to  match  the  peak  of 
the  brightest  point  In  the  Image,  and  the  brightness  and  position  of  the 
peak  was  noted.  The  second  and  third  points  were  CLEANed  in  a  similar 
fashion  In  succession.  Figure  7-l(d)  shows  the  brightness  and  positions 
of  the  three  CLEANed  peaks.  The  brightness  and  positions  are  exactly 
the  same  (to  within  an  overall  translation)  as  the  brightness  and 


Figure  7- 


Shift-and-Add  Algorithm  for  a  Point  Object,  (a)  Original 
object;  (b)  image  blurred  by  atmospheric  turbulence; 

(c)  output  image  from  shi ft-and-add;  (d)  CLEANed  version 
of  output  image. 
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positions  of  the  points  In  the  original  object.  Although  no  noise  Mas 
present  In  the  data  (but  only  156  blurred  Images  were  used),  this  Is  a 
very  Impressive  result,  demonstrating  the  power  of  the  shlft-and-add 
method  for  this  type  of  object. 

Figure  7-2  shows  a  similar  experiment  for  an  extended  object.  The 
object,  shown  in  Figure  7-2(a),  was  purposely  chosen  to  be  one  that 
satisfies  the  requirement  that  it  have  bright  delta-function-like 
components.  That  Is,  It  should  be  one  of  the  easier  extended  objects 
for  shlft-and-add  to  reconstruct.  Two  examples  of  the  blurred  Images  of 
the  object  are  shown  In  Figures  7-2(b)  and  (c).  The  result  of  using 
shlft-and-add  on  156  blurred  Images  is  shown  In  Figure  7-2(d).  In  this 
case  there  Is  very  little  Information  about  the  original  object  In  the 
shlft-and-add  Image.  Further  processing  using  the  CLEAN  method  did  not 
Improve  the  result. 

This  one  set  of  experiments  was  not  sufficient  to  fully  delineate 
the  types  of  objects  for  which  shlft-and-add  Is  effective,  but  we  did 
demonstrate  that  shlft-and-add  works  very  well  for  an  object  consisting 
of  a  small  number  of  delta-functlon-llke  points  dominated  by  a  brightest 
one,  but  works  very  poorly  for  an  extended  object,  even  one  containing 
Isolated  bright  points. 
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RECONSTRUCTION  OF  COMPLEX-VALUED  OBJECTS  USING  SUPPORT  CONSTRAINT 


As  discussed  elsewhere  In  this  report,  for  the  astronomy  problem 
one  has  In  the  object  domain  a  nonnegativity  constraint  and  a  weak 
(loose)  support  constraint.  There  are  however  some  problems  for  which 
the  object  Is  complex-valued  or  bipolar,  precluding  a  nonnegativity 
constraint,  but  for  which  the  support  constraint  Is  much  stronger 
(tighter).  The  same  Iterative  reconstruction  algorithm  as  for  the 
astronomy  problem  can  be  used,  only  without  applying  the  nonnegativity 
constraint.  From  the  theory  of  Bruck  and  Sodin  [27],  one  might  expect 
the  solution  to  usually  be  unique.  For  supports  known  to  have  certain 
shapes,  such  as  a  triangular  shape,  the  solution  Is  known  to  be  unique 
(see  Appendix  B). 

Reconstruction  experiments  using  only  a  support  constraint  were 
performed  on  objects  having  various  support  constraints  to  test  the 
Importance  of  different  types  of  support  constraints.  Figure  8-l(a) 
shows  an  object  having  triangular  support  and  nonzero  values  in  Its 
three  corners.  These  conditions  ensure  that  the  object  is  uniquely 
related  to  Its  Fourier  modulus  (see  Appendix  B).  The  Image  shown  In 
Figure  8-l(b)  was  reconstructed  from  the  Fourier  modulus  and  the  a 
priori  knowledge  of  the  triangular  support  using  the  iterative  Fourier 
transform  algorithm.  Nonnegativity  was  not  used  although  the  object 
happens  to  be  nonnegative.  In  this  case  the  support  was  known  very 
precisely.  The  algorithm  converged  very  rapidly  to  the  correct 
solution. 

Since  the  uniqueness  proof  requires  the  three  corners  to  be 
nonzero,  we  wanted  to  determine  the  Importance  of  nonzero  corners  to  the 
Iterative  algorithm.  The  same  experiment  was  performed  for  the  object 
shown  In  Figure  8-l(c),  which  Is  Identical  to  the  object  shown  In  Figure 
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Figure  8-1.  Examples  of  Reconstruction  from  Fourier  Modulus. 

(a)  Object  in  triangle  with  bright  corners,  (b)  recon¬ 
structed  image;  (c)  object  in  triangle  with  zeroed  corn¬ 
ers,  (d)  reconstructed  image;  (e)  object  in  triangle  with 
tapered  edges,  (f)  reconstructed  image. 
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8- 1 ( a )  but  with  the  corners  zeroed  out.  The  Image  shown  in  Figure 
8-l(d)  was  reconstructed  from  its  Fourier  modulus  and  Is  the  correct 
solution,  but  more  iterations  were  required  In  this  case  than  for  the 
object  having  three  bright  corners.  Therefore  the  brightness  of  the 
corners  does  play  a  role,  but  not  a  crucial  one.  The  effect  of  the 
sharpness  of  the  edges  of  the  object  was  also  Investigated.  A  third 
object,  shown  In  Figure  8-l(e),  which  is  identical  to  the  object  shown 
In  Figure  8- 1  ( c )  except  that  Its  edges  are  tapered  rather  than  being 
abrupt,  was  formed.  The  Image  resulting  after  over  a  hundred  iterations 
of  the  Iterative  Fourier  transform  algorithm  is  shown  in  Figure  8- 1(f) . 
Although  the  Image  is  very  recognlzble.  It  has  a  noisy  appearance.  The 
iterative  algorithm  found  it  much  more  difficult  to  reconstruct  this 
Image  than  the  ones  with  abrupt  or  sharp  edges.  Therefore  it  appears 
that  edges  (although  not  absolutely  essential)  are  very  Important  to  the 
ability  of  the  Iterative  algorithm  to  reconstruct  images  using  only  a 
support  constraint  In  the  object  domain. 

Under  an  Internally  funded  ERIM  program,  these  results  were 
extended  to  the  case  of  complex- valued  objects  for  applications  such  as 
synthetic-aperture  radar  (SAR)  [25].  Since  those  results  are  pertinent 
here,  they  will  be  briefly  reviewed.  The  Idea  Is  to  have  a  SAR  sensor 
that  does  not  require  an  accurate  local  oscillator,  phase-coherent  chirp 
pulses,  or  compensation  of  sensor  platform  motion.  This  could  be  done 
If  one  could  reconstruct  the  Image  without  accurate  knowledge  of  the 
phase  of  the  SAR  signal  history.  This  might  be  possible  by  using  the 
Iterative  Fourier  transform  algorithm  If  a  strong  support  constraint 
were  present.  One  might  have,  for  example,  the  ability  to  Illuminate 
the  target  area  with  an  lllumlnaton  pattern  of  known  shape  or  have  the 
far-fleld  pattern  of  the  receive  antenna  accept  reflected  radiation  from 
an  area  of  known  shape.  The  modulus  Information  together  with  the 
support  (beam  shape)  constraint  could  then  be  combined  by  the  iterative 
Fourier  transform  algorithm  to  arrive  at  a  diffraction-limited  Image. 


Figure  8-2  shows  an  example  of  a  reconstruction  experiment  [25]  of 
this  type.  Figure  8- 2(a)  shows  the  magnitude  of  a  64  x  64  pixel 
sub-area  of  a  complex-valued  SEASAT  SAR  image  of  an  area  of  land.  A 
binary  mask  was  formed  to  define  the  Illumination  pattern  of  a 
hypothetical  antenna.  The  Illumination  pattern  consists  of  a  pair  of 
ellipses,  each  of  3:1  aspect  ratio.  A  pattern  consisting  of  two 
separated  parts  was  chosen  because  theory  indicates  that  the  solution  is 
more  likely  to  be  unique  In  that  case  [13].  The  object,  the  magnitude 
of  which  Is  shown  In  Figure  8-2(b),  was  obtained  by  taking  the  product 
of  the  Image  shown  In  Figure  8-2(a)  and  the  illumination  pattern.  The 
modulus  of  the  Fourier  transform  of  the  object  was  computed  and  is  shown 
in  Figure  8-2(c).  This  is  equivalent  to  the  modulus  of  the  SAR  signal 
history  that  would  have  been  collected  had  the  terrain  been  Illuminated 
by  the  fixed  illumination  pattern  consisting  of  the  two  ellipses.  The 
Fourier  modulus  was  then  used  together  with  the  known  illumination 
pattern  as  a  support  constraint  to  reconstruct  the  image  which  is  shown 
In  Figure  8-2(d).  The  reconstruction  was  essentially  perfect. 

Figure  8-3  shows  further  examples  of  similar  reconstruction 
experiments  performed  under  the  current  effort.  The  goal  of  this  set  of 
experiments  was  to  explore  the  effects  of  employing  different  types  of 
support  constraints.  Specifically,  we  wanted  to  explore  what  effect  the 
separation  of  the  parts  of  the  support  had  on  the  success  of  the 
Iterative  reconstruction  algorithm.  As  shown  in  Figure  8-2,  for  widely 
separated  support  parts,  the  iterative  algorithm  performed  very  well. 
Figure  8-3(a)  and  8-3(b)  show  the  same  object  and  reconstructed  image, 
respectively,  but  at  a  different  scale.  Figures  8-3(c)  and  8-3(d)  show 
a  second  object  and  reconstructed  Image,  respectively,  for  a  case  in 
which  the  separation  between  the  two  parts  of  the  support  is  much 
smaller.  Again  the  reconstructed  image  is  very  faithful.  Figure  8-3(e) 


Figure  8-2.  Example  of  Reconstructing  a  Complex-Valued  SAR  Image  from 
the  Modulus  of  Its  Fourier  Transform  Using  an  Illumina¬ 
tion-Pattern  Support  Constraint  (a  Pair  of  Ellipses). 

(a)  Magnitude  of  terrain  image  with  broad  illumination 
pattern;  (b)  magnitude  of  ideal  terrain  image  with  special 
illumination  pattern;  (c)  Fourier  (phase  history)  modulus; 
(d)  magnitude  of  reconstructed  image.  (Taken  from  [ 25 ] . ) 
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Figure  8-3.  Examples  of  Reconstructing  a  SAR  Image  from  Modulus  of  Its 
Fourier  Transform  Using  Various  Support  Constraints  with 
the  Iterative  Transform  Algorithm.  Illuminated  objects: 
(a),  (c),  (o),  (g);  respective  reconstructed  images:  (b), 
(d),  (f),  (;.). 
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and  8-3(f)  shows  a  third  object  and  reconstructed  image,  respectively, 
for  a  case  in  which  the  two  parts  overlap  (that  is,  the  support  is 
contiguous).  In  this  case  the  reconstructed  image,  even  after  several 
hundred  iterations,  does  not  closely  resemble  the  object,  although  upon 
close  inspection  one  can  find  some  features  in  common.  For  a  fourth 
case  (not  shown)  in  which  the  support  consisted  of  a  single  ellipse,  the 
iterative  reconstruction  algorithm  did  not  produce  a  recognizable  image 
after  several  hundred  iterations.  Figure  8-3(g)  and  8-3(h)  show  a  fifth 
object  and  reconstructed  image,  respectively,  for  a  case  in  which  the 
support  was  shaped  like  a  donut  having  a  hole  offset  from  the  center. 
In  this  case  the  reconstructed  image  is  very  faithful.  This  last 
example  is  curious  since  it  does  not  have  separated  parts,  yet  a 
one-dimensional  cut  through  the  center  of  the  support  does  have 
separated  parts. 

The  examples  of  Figure  8-3  demonstrate  that  the  separated  nature  of 
the  support  does  have  an  important  effect  on  the  succes  of  the  iterative 
reconstruction  algorithm,  and  that  further  investigations  along  these 
lines  is  warranted. 

The  ability  demonstrated  above  to  reconstruct  a  complex-valued 
image  from  the  modulus  of  its  Fourier  transform  using  only  a  support 
constraint  may  have  very  Important  implicatons  for  fine-resolution 
coherent  imaging  systems  such  as  SAR. 
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OF  AN  ITERATIVE  ALGORITHM 
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Transformations  in  Optical  Signal  Processing,  Proc.  SPIE  373,  147-160 
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Reconstruction  and  synthesis  applications  of  an  iterative 
algorithm 


J.  R.  Fienup  Abstract.  This  paper  reviews  the  Gerchberg-Saxton  algorithm  and  varia- 

Environmental  Research  Institute  of  tions  thereof  that  have  been  used  to  solve  a  number  of  difficult  recon- 

Michigan  struction  and  synthesis  problems  in  optics  and  related  fields.  It  can  be 

P  O  Box  8618  used  on  any  problem  in  which  only  partial  information  (including  both 

Ann  Arbor,  Michigan  48107  measurements  and  constraints)  of  the  wavefront  or  signal  is  available  in 

one  domain  and  other  partial  information  is  available  in  another  domain 
(usually  the  Fourier  domain).  The  algorithm  combines  the  information  in 
both  domains  to  arrive  at  the  complete  description  of  the  wavefront  or 
signal.  Various  applications  are  reviewed,  including  synthesis  of  Fourier 
transform  pairs  having  desirable  properties  as  well  as  reconstruction 
problems.  Variations  of  the  algorithm  and  the  convergence  properties  of 
the  algorithm  are  discussed. 


1.  INTRODUCTION 

There  exist  many  problems  that  are  very  difficult  to  solve  in 
asironomy,  x-ray  crystallography,  electron  microscopy,  spec¬ 
troscopy,  wavefront  sensing,  holography,  particle  scattering, 
superresolution,  radar  signal  and  antenna  synthesis,  filter  design, 
and  other  disciplines  that  share  an  important  feature.  These  are 
problems  that  involve  the  reconstruction  or  synthesis  of  a 
wavefront  (or  an  object  or  a  signal,  etc.)  when  partial  information 
or  constraints  exists  in  each  of  two  different  domains.  The  second 
domain  is  usually  the  Fourier  transform  domain.  This  paper 
describes  a  method  of  combining  all  the  available  information  in 
the  two  domains  to  arrive  at  a  complete  description,  thereby  solv¬ 
ing  the  problems. 

The  problems  fall  into  two  general  categories:  (1)  reconstruct  the 
entire  information  about  a  function  (an  image,  wavefront,  signal, 
etc.)  when  only  partial  information  is  available  in  each  of  two  do¬ 
mains;  and  (2)  synthesize  a  (Fourier)  transform  pair  having 
desirable  properties  in  both  domains.  A  reconstruction  problem 
arises  when  only  partial  information  is  measured  in  one  domain, 
and  in  the  other  domain  either  partial  information  is  measured  or 
certain  constraints  are  know  n  a  prion.  The  information  available  in 
any  one  domain  is  insufficient  to  reconstruct  the  function  or  its 
transform.  A  synthesis  problem  typically  arises  when  one  wants  the 
transform  of  a  function  to  have  certain  desirable  properties  (such  as 
uniform  spectrum,  low  sidelobes,  etc.)  while  the  function  itself 
must  satisfy  certain  constraints  or  have  certain  desirable  properties. 
Because  arbitrary  sets  of  properties  and  constraints  can  be  con¬ 
tradictory,  there  may  not  exist  a  transform  pair  that  is  completely 
desirable  and  satisfies  all  the  constraints.  Nevertheless,  one  seeks  a 
transform  pair  that  comes  as  close  as  possible  to  having  the 
desirable  properties  and  satisfy  ing  the  constraints  in  both  domains. 

Both  the  reconstruction  and  the  synthesis  problems  can  be  ex¬ 
pressed  as  follows,  if  the  meaning  of  the  word  “constraints-'  is 
broadened  to  include  any  kind  of  measured  data,  desirable  proper¬ 


ties,  or  a  priori  conditions: 

(divert  a  set  of  constraints  placed  on  a  function  and  another 
set  of  constraints  placed  on  its  transform,  find  a  transform 
pair  (i.e.,  a  function  and  its  transform)  that  satisfies  both  sets 
of  constraints. 

Once  a  solution  is  found  to  such  a  problem,  the  question  often 
remains:  is  the  solution  unique?  For  synthesis  problems,  the 
uniqueness  is  usually  unimportant— one  is  satisfied  with  any  solu¬ 
tion  that  satisfies  all  the  constraints:  often  a  more  important  prob¬ 
lem  is  whether  there  exists  any  solution  that  satisfies  what  may  be 
arbitrary  and  conflicting  constraints.  For  reconstruction  problems, 
the  uniqueness  properties  of  the  solution  are  of  central  importance. 
If  many  different  functions  satisfy  ing  the  constraints  could  give  rise 
to  the  same  measured  data,  then  a  solution  that  is  found  could  not 
be  guaranteed  to  be  the  correct  solution.  The  question  of  unique¬ 
ness  must  be  studied  for  each  problem.  Fortunately,  as  will  be 
described  later,  for  some  important  reconstruction  problems  the 
solution  usually  is  unique. 

An  effective  approach  to  solving  the  large  class  of  problems 
described  above  is  the  use  of  iterative  algorithms  related  to  the 
Gerchberg-Saxton  algorithm.1  The  algorithms  involve  the  iterative 
transformation  back  and  forth  between  the  two  domains,  with  the 
known  constraints  applied  repetitively  in  each  domain. 

The  basic  algorithm  is  presented  in  Sec.  2.  A  number  of  different 
applications  having  different  types  of  constraints  are  described, 
and  examples  are  shown  in  Sec.  In  See.  4  the  convergence  prop¬ 
erties  of  the  algorithm  are  discussed,  and  improved  versions  of  the 
algorithm  arc  reviewed.  A  brief  summary  and  comments  arc  in¬ 
cluded  in  Sec.  5. 

2.  THK  BASIC  ITKRATIVK  ALGORITHM 

The  first  published  account  of  the  iterative  algorithm  was  iis  use  bv 
(ierchberg  and  Saxton1  to  solve  the  electron  microscopy  problem. 
For  this  problem  both  the  modulus  (magnitude)  of  a  complex- 
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valued  image  and  the  modulus  of  its  Fourier  transform  are 
measured,  and  the  goal  is  to  reconstruct  the  phase  in  both  domains. 
Apparently  unknown  to  Gerchberg  and  Saxton,  the  method  was  in¬ 
vented  somewhat  earlier  by  Hirsch,  Jordan,  and  l.esem2  to  solve  a 
synthesis  problem  for  computer-generated  holograms  that  has  a 
similar  set  of  constraints.  (This  will  be  described  later  in  more 
detail.)  The  method  was  again  reinvented  for  a  similar  problem  in 
computer  holography  by  Gallagher  and  Liu.-'  The  fact  that  the 
algorithm  was  invented  repeatedly  testifies  to  its  simplicity  and  ef¬ 
fectiveness. 

2.1.  Gerchberg-Saxton  algorithm 

In  what  immediately  follows,  the  iterative  algorithm  is  described  in 
terms  of  its  application  to  the  electron  microscopy  reconstruction 
problem.  An  excellent  treatment  of  the  electron  microscopy  phase 
problem  and  its  solution  by  this  and  other  methods  can  be  found  in 
Ref.  4.  Later  it  is  shown  how  to  apply  the  same  principles  to  a  large 
class  of  problems. 

Suppose  that  the  electron  wave  function  in  an  image  plane  is 
described  by  the  two-dimensional  (2-D)  complex-valued  function 
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Fig.  1.  Block  diagram  of  the  iterative  error-reduction  algorithm. 


f(x)  —  f(x)  c,Clx)  .  (I) 

Its  Fourier  transform,  the  wave  function  in  a  far-field  diffraction 
plane,  is  given  by 

OO 

F(u)  =  F(u)  eiS(u>  =  #[f(x)|  =  f  f(x)ei’’ru”ldx,  (2) 

-  OO 

where  x  and  u  are  the  vector  coordinates  in  the  spatial  (image)  do¬ 
main  and  the  spatial  frequency  (far-field  diffraction)  domain, 
respectively.  The  notation  used  throughout  this  paper  is  that  func¬ 
tions  represented  by  capital  letters  are  the  Fourier  transforms  of  the 
functions  represented  by  the  corresponding  lower-case  letters.  It  is 
assumed  that  the  intensity  spatial  distributions  are  measured  in 
each  domain,  but  the  phase  information  is  lost.  Therefore,  one 
wishes  to  reconstruct  \i(x)  and  fl(x)  from  f(x)  and  F(u)  . 

The  iterative  algorithm  for  solving  this  problem  is  depicted  in 
Fig.  1 .  One  iteration  (the  k,h  iteration)  of  the  algorithm  proceeds  as 
follows.  A  trial  solution  tor  the  wave  function  (an  estimate  of  the 
wave  function),  gk(x),  is  Fourier  transformed  yielding 

Gk(u)  -  Gk(u)  exp[iok(u)|  ,<*  lgk(x)|  .  (3) 

Then  a  new  Fourier-domain  function,  Gk(u),  is  formed  by  replac¬ 
ing  the  computed  Fourier  modulus  by  the  measured  Fourier 
modulus,  F(u)  ,  and  keeping  the  computed  phase: 

Gk(u)  -  Flu)  exp[iok(u)|  (4) 

The  resulting  Gk(u),  which  is  in  agreement  with  all  the  known 
measurements  and  constraints  in  the  Fourier  domain,  is  inverse 
Fourier  transformed,  yielding  the  wave  function  gk(x).  The  itera¬ 
tion  is  completed  by  forming  a  new  estimate  for  the  wave  function, 
gk  .  |(x),  which  is  obtained  by  replacing  the  computed  modulus  of 
gk(x)  with  the  measured  modulus  f(x)  ,  and  keeping  the  computed 
phase. 

The  algorithm  consists  of  no  more  than  enforcing  what  informa¬ 
tion  is  available  on  the  wave  function,  Fourier  transforming,  im¬ 
posing  what  information  is  available  on  the  wave  function’s 
Fourier  transform,  inverse  transforming,  and  repeating  these  sim¬ 
ple  operations  for  a  number  of  iterations.  What  makes  the 
algorithm  practical  is  the  existence  of  a  fast  Fourier  transform' 
(FFT),  so  that  the  number  of  computations  per  iteration  goes  only 
as  NIogN,  where  N  is  the  number  of  samples  of  the  function  com¬ 
puted.  This  compares  very  favorably  with  some  other  iterative 


methods,  such  as  Newton-Raphson,4  for  which  the  number  of  com¬ 
putations  per  iteration  goes  as  N\ 

A  measure  of  the  progress  of  the  iterations,  and  a  criterion  by 
which  one  can  determine  when  a  solution  has  been  found,  is  the 
normalized  mean-squared  error,  which  is  defined  in  the  Fourier  do¬ 
main  by  _ 

OC 

J  [  Gk(u)  -  Flu)  ]■  du 

¥-■»  -OO 

Ef  = -  (5) 

OO 


-OO 


or  in  the  image  domain  by 


Eo 


00 

J  I  gg(x)  -  f(x)  ):dx 
-00 


00 

j  f(x)  -dx 

-  OO 


(6) 


It  has  been  shown  that  the  algorithm  converges  in  the  sense  that  the 
mean-squared  error  can  only  decrease  at  each  iteration.1  The 
issue  of  convergence  will  be  discussed  in  greater  detail  m  Sec.  4. 


2.2.  Krror-reduclion  iterative  algorithm 

It  is  now  known  that  with  slight  modifications  this  same  algorithm 
can  be  applied  to  many  different  problems  having  a  variety  of 
available  constraints  or  measurements.'  Let  the  function  f(x)  repre¬ 
sent  a  wavefront,  an  object,  a  signal,  an  antenna  array,  a  spectral 
density  function,  an  electron  density  function,  etc.,  where  \  is  an 
N-dimensional  vector  (spatial,  angular,  time,  etc.)  coordinate. 
Depending  on  the  problem,  f(x)  may  be  complex  valued  or  real 
valued  and,  if  real,  may  or  may  not  be  nonnegative  Its  Fourier 
transform,  F(u),  is  given  by  Eq.  (2)  and  is  complex  valued  for  most 
problems.  The  N-dimensional  vector  u  is  a  (spatial,  angular,  time, 
etc.)  frequency  coordinate.  One  can  instead  consider  another 
transformation  of  f(x),  such  as  the  Fresnel  transform,  which  has 
been  used  for  more  than  one  problem. ;  s  u  l  or  simplicity  of  discus¬ 
sion,  the  Fourier  transform  will  be  assumed,  but  the  reader  should 
keep  in  mind  that  what  is  said  also  applies  to  a  number  of  other 
transformations  as  well  (although  the  method  becomes  less  attrac- 
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live  if  a  fast  transform  algorithm  is  not  available). 

With  only  slight  modifications,  the  Gerchberg-Saxton  algorithm 
can  be  used  to  solve  the  wide  class  of  problems  described  in  Sec.  1 . 
Referring  again  to  the  block  diagram  of  the  algorithm  in  Fig.  1,  all 
that  is  required  is  to  impose  constraints  in  each  domain  that  are 
pertinent  to  the  problem  of  interest  At  the  k,h  iteration,  gk(x),  an 
estimate  of  f(x),  is  Fourier  transformed,  yielding  Gk(u),  which  is 
given  by  Eq.  (3).  Then  a  new  Fourier-domain  function  Gk(u)  is 
formed  from  Gk(u)  by  making  the  smallest  possible  changes  in 
Gk(u)  that  allow  it  to  satisfy  the  Fourier-domain  constraints.  For 
example,  if  the  Fourier-domain  constraint  is  that  the  Fourier 
modulus  equals  F(u)  over  some  region  of  the  Fourier  domain, 
then  F(u)  is  substituted  for  i  Gk(u)  |  in  that  region.  The  new 
Fourier-domain  function  Gk(u),  which  satisfies  the  Fourier-domain 
constraints,  is  inverse  Fourier  transformed  to  yield  gk(x).  To  com¬ 
plete  one  iteration,  a  new  estimate  gk  +  ^x)  is  formed  from  gk(x)  by 
making  the  smallest  possible  changes  in  gk(x)  that  allow  it  to  satisfy 
the  function-domain  constraints.  One  example  is  that  if  the  func¬ 
tion  is  complex  valued  and  it  is  constrained  to  have  a  modulus 
equal  to  f(x)  over  some  region  of  space,  then  :  f(x)|  is  substituted 
for  gk(x)  in  that  region.  A  special  case  of  this  is  when  the  func¬ 
tion  is  to  be  zero  outside  a  certain  interval  (the  Fourier  function  is 
bandlimited).  Another  example  is  that  if  the  function  is  constrained 
to  be  nonnegative,  then  gk  +  ,(x)  is  set  equal  to  gk(x)  for  those  x 
where  gk(x)  >  0,  and  gk  *  ](x)  is  set  equal  to  zero  for  those  x  where 
g.'  fx)  <  0.  In  summary,  one  transforms  back  and  forth  between  the 
two  domains,  forcing  the  function  to  satisfy  the  constraints  in  each 
domain. 

For  reconstruction  problems,  whatever  characteristics  of  the  ac¬ 
tual  F(u)  and  f(x)  that  are  measured  or  are  known  a  priori  are  im¬ 
posed  on  Gk(u)  and  gk(x),  respectively.  For  synthesis  problems, 
one  imposes  on  Gk(u)  and  gk(x)  whatever  characteristics  one  might 
desire  F(u)  and  f(x),  respectively,  to  have.  Once  the  constraints  are 
defined,  the  algorithm  proceeds  the  same  for  synthesis  problems  as 
for  reconstruction  problems.  In  fact,  there  are  some  synthesis  prob¬ 
lems  that  are  mathematically  indistinguishable  from  some 
reconstruction  problems,  and  they  are  handled  identically  by  the 
algorithm. 

The  first  iteration  of  the  algorithm  can  be  started  in  a  number  of 
ways,  for  example,  by  setting  g,(x)  or  <6,(x)  equal  to  an  array  of 
random  numbers.  The  iterations  continue  until  a  Fourier  transform 
pair  is  found  that  satisfies  all  the  constraints  in  both  domains  to 
within  the  desired  accuracy  (or,  if  convergence  is  too  slow,  until 
one  loses  interest  or  the  money  runs  out).  The  mean-squared  error 
can  generally  be  defined  in  the  Fourier  domain  by 
00 

J  Gk(u)  -  Gk(u)  :du 

F.f  - - --  P) 

OO 

Gk(u)  :  du 

-  OO 

or  in  the  function  domain  by 

OC 

f  gk  .  |(G  -  gk(x)  :dx 

!•  *  <«> 

OC 

j  gk(x)  ;  d\ 

-  OC 

In  each  of  these  two  expressions,  the  integrand  in  the  numerator  is 
the  squared  modulus  ol  the  amount  by  which  the  computed  func¬ 
tion  violates  the  constraints  in  that  domain.  It  is  easily  seen  that 


these  expressions  reduce  to  Eqs.  (5)  and  (6),  respectively,  for  the 
electron  microscopy  problem. 

Just  as  in  the  electron  microscopy  problem,  for  problems  having 
other  sets  of  constraints  it  will  be  shown  in  Sec.  4  that  the  algorithm 
converges,  that  is,  the  error  decreases  at  each  successive  iteration. 
The  algorithm  depicted  in  Fig,  1  may  be  referred  to  as  the  “error- 
reduction”  algorithm  for  that  reason,  as  well  as  to  distinguish  it 
from  algorithms  described  in  Sec.  4  that  are  related  to  it  but  con¬ 
verge  faster.  Typically,  the  error  is  reduced  very  rapidly  for  the  first 
few  iterations  of  the  error-reduction  algorithm,  but  more  slowly  for 
later  iterations.  For  some  applications,  the  error-reduction 
algorithm  has  been  very  successful  in  finding  solutions  using  a 
reasonable  number  of  iterations.  However,  for  some  other  applica¬ 
tions,  the  mean-squared  error  decreases  extremely  slowly  with  each 
iteration,  and  an  impractically  large  number  of  iterations  is  re¬ 
quired.  The  improved  algorithms  described  in  Sec.  4  do  much  to 
alleviate  this  problem. 

2.3.  Alternative  descriptions  of  the  algorithm 

Once  a  solution  (i.e.,  a  Fourier  transform  pair  satisfying  all  the 
constraints  in  both  domains)  is  found,  the  error-reduction 
algorithm  ceases  to  make  changes  to  the  estimate,  and  the 
algorithm  locks  on  to  the  solution.  The  operations  of  enforcing  the 
constraints  in  each  domain  would  then  leave  the  function  estimate 
and  its  Fourier  transform  unaltered,  since  they  already  satisfy  the 
constraints.  Now  let  u.s  define  the  operation  Sig(\)J  as  the  suc¬ 
cessive  Fourier  transformation  of  g(.x),  followed  by  the  imposition 
of  the  Fourier  domain  constraints,  followed  by  inverse  Fourier 
transformation,  followed  by  imposition  of  the  object  domain  con¬ 
straints.  That  is,  the  operation  S  is  just  the  performance  of  one 
iteration  of  the  error-reduction  algorithm,  and 

gk  t  ,(x)  =  S[gk(x)]  .  (9) 

From  the  discussion  above,  it  is  evident  that  any  solution  f(\)  must 
satisfy  the  relation 

f(x)  =  S[f(x)] .  GO) 

When  presented  in  this  form,  it  is  seen  that  the  error-reduction 
algorithm  is  a  particular  implementation  of  the  method  of  suc¬ 
cessive  approximations.11’ 

The  method  of  successive  approximations  can  be  more  easily 
understood  from  the  following  simple  example.  Suppose  one 
wishes  to  solve  the  following  equation  for  y: 

4yJ  -  4y  *  I  -  0  .  0  1 1 

Based  on  the  relation  y  -  yJ  4  1  4,  one  could  write 

yk  •  i  =  si«>k>  *  1  4  <l2> 

Using  the  method  of  successive  approximations  to  find  the  solu¬ 
tion,  one  would  pick  an  initial  estimate,  say  v0  -■  0.1,  and  employ¬ 
ing  F.q.  (12)  compute  y,  -  0.2501,  y;  -  0.2539,  etc.,  and  rapidly 
converge  to  the  solution  y  0.2541737  ....  However,  it  con¬ 
verges  to  y  only  for  y(,  <  y  ”  -  0.8 967902  ....  For  y0  >  y ",  Eq. 
(12)  diverges;  and  for  y0  y  ",  it  stays  at  y ",  the  second  solution. 
On  the  other  hand,  one  could  just  as  logically  have  chosen 

yk  .  ,  S:(yk)  (yk  -  I  4)1  4  .  (13) 

This  second  form  converges  to  the  second  solution  y  for  y0  >  y  . 
diverges  for  y0  <  v  ,  and  stays  at  y  for  y0  v  .  Figure  2,  a 
graphical  representation  of  Eq.  (12),  shows  the  two  solutions,  v 
and  y".  The  irregular  staircase  between  the  two  curves  v  and  >J  • 
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Fig.  2.  Method  of  successive  approximations  for  solving  4y4  -  4y 
+  1=0. 


1/4  indicates  how  the  estimate  yk  approaches  the  two  solutions. 
Criteria  on  the  derivative  of  S(y)  determine  whether  the  algorithm 
converges." 

The  error-reduction  algorithm,  as  described  by  Eqs.  (9)  and  (10), 
is  analogous  to  the  example  of  successive  approximations  described 
above,  except  that  instead  of  operating  on  a  scalar  y,  it  operates  on 
a  function  g(x).  As  seen  from  the  example,  the  method  of  suc¬ 
cessive  approximations  may  or  may  not  converge,  depending  on  the 
particular  form  chosen  and  on  the  initial  estimate.  Fortunately,  as 
will  be  discussed  further  in  Sec.  4,  the  error-reduction  algorithm 
never  diverges.  It  may,  however,  stagnate.  A  simple  example  of 
stagnation  of  the  method  of  successive  approximations  is  shown  by 
the  following.  In  solving  x  =  2  -  x  (which  has  the  obvious  solution 
x  =  1),  starting  with  the  initial  estimate  x0,  one  obtains  x,  =  2-x0, 

x2  =  2  -  (2  -  x0)  =  x0 . x2k.,  =  2  -  x0,  x2k  =  x0,  etc.,  and  no 

progress  is  made  toward  the  solution. 

Another  way  of  understanding  the  error-reduction  algorithm, 
applicable  for  certain  sets  of  constraints,  is  the  alternating  projec¬ 
tion  of  the  function  onto  specified  subspaces  in  a  Hilbert  space.12 
This,  along  with  the  possibility  of  closed-form  solutions,13  is 
discussed  in  the  contribution  to  this  volume  by  Marks  and  Smith. 

3.  APPLICATIONS 

A  large  number  of  important  problems  in  optics  and  related  fields 
fit  the  problem  description  in  Sec.  1  and  can  be  solved  by  the 
iterative  algorithm  (by  the  error-reduction  algorithm  described  in 
Sec.  2  and  the  related  algorithms  described  in  Sec.  4).  One  par¬ 
ticular  application,  that  of  spectral  extrapolation  or  superrcsolu- 
tion,  is  discussed  in  detail  in  the  contribution  to  this  volume  by 
Marks  and  Smith.  In  this  section,  several  classes  of  applications  are 
listed,  followed  by  more  detailed  discussions  of  some  of  the  ap¬ 
plications,  including  examples. 

In  Sec.  I,  a  distinction  was  made  between  reconstruction  prob¬ 
lems  and  synthesis  problems.  Another  useful  way  to  classify  such 
problems  is  according  to  the  type  of  information  available.  For  one 
set  of  problems,  the  modulus  (magnitude  or  amplitude)  of  a 
complex-valued  function  and  the  modulus  of  its  Fourier  transform 


are  measured  (or  are  given),  and  one  wishes  to  know  the  phase  of 
the  Fourier  transform  pair  in  both  domains.  These  include  the 
phase  retrieval  problem  in  electron  microscopy,  the  phase  retrieval 
problem  in  wavefront  sensing,  the  design  optimization  of  radar 
signals  and  antenna  arrays  having  desirable  properties,  and  phase 
coding  and  spectrum  shaping  problems  for  computer-generated 
holograms  and  other  applications.  These  applications  often  involve 
the  Fresnel  transform  for  the  near-field  case  instead  of  the  Fourier 
transform. 

For  another  set  of  problems,  the  function  is  known  to  be  real  and 
nonnegative  and  the  modulus  of  its  Fourier  transform  is  measured. 
These  include  the  phase  problems  of  x-ray  crystallography,  Fourier 
transform  spectroscopy,  imaging  through  atmospheric  turbulence 
using  interferometer  data,  and  pupil  function  determination. 

For  another  set  of  problems,  a  low-resolution  (i.e.,  a  low-pass 
filtered)  version  of  a  function  is  measured  (i.e.,  its  complex  Fourier 
transform  is  measured  only  over  a  certain  interval),  and  the  func¬ 
tion  is  known  to  have  a  finite  extent  (i.e.,  it  is  zero  outside  of  some 
known  region  of  support).  This  is  the  spectral  extrapolation  or 
superresolution  problem  for  band-limited  time  signals  or  for  im¬ 
aging  of  objects  of  finite  extent. 

For  another  set  of  problems,  the  function  is  known  to  be  non¬ 
negative  and  of  finite  extent  and  its  complex  Fourier  transform  is 
measured  only  over  a  partially  filled  aperture.  These  include  the  in¬ 
terpolation  of  the  complex  visibility  function  for  long  baseline 
radio  interferometry  and  the  missing-cone  problem  in  x-ray 
tomography. 

For  still  another  set  of  problems,  the  modulus  of  a  complex¬ 
valued  function  is  given,  and  one  wishes  to  find  an  associated  phase 
function  that  results  in  a  Fourier  transform  whose  complex  values 
fall  on  a  prescribed  set  of  quantized  complex  values.  These  include 
the  reduction  of  quantization  noise  in  computer-generated 
holograms  and  in  coded  signal  transmission. 

Another  problem  is  to  reconstruct  the  modulus  of  a  complex¬ 
valued  function  from  the  phase  of  the  function,  given  the  fact  that 
the  Fourier  transform  of  the  function  has  finite  support. 

The  number  of  types  of  problems  solvable  by  the  iterative 
algorithm  appears  to  be  limited  only  by  one’s  ingenuity  in  defining 
different  combinations  of  information  that  might  be  available  in 
each  of  two  domains. 

3.1.  Modulus— modulus  constraints 
3.  /.  /.  Electron  microscopy 

Among  the  applications  for  which  the  modulus  is  given  in  each  of 
two  domains,  the  electron  microscopy  phase  retrieval  problem  was 
one  of  the  earliest  applications  of  the  error-reduction  algorithm 
and  has  been  the  problem  most  heavily  investigated. '-4*-14"  The 
error-reduction  (Gerchberg-Saxton)  algorithm  has  been  shown  to 
perform  very  successfully  for  this  problem,  and  the  solution  is 
usually  unique."  The  reader  is  referred  to  a  book  by  Saxton4  for  a 
thorough  review. 

3.1.2.  Spectrum  shaping 

A  second  application  for  which  the  modulus  is  given  in  each  of  two 
domains  is  the  spectrum  shaping  problem.  Spectrum  shaping  is  a 
synthesis  problem  that  can  be  stated  as  follows:  given  the  modulus 
f(x)  of  a  complex-valued  wavefront,  g(x)  =  f(x)  exp[ift(x)], 
find  a  phase  function  0(x)  such  that  .^[g(x)|  is  equal  to  a  given 
spectrum  F(u)  .  Such  a  problem  is  the  one  suggested  by  the  Escher 
engraving  shown  in  Fig.  3,  in  which  a  bird  transforms  into  a  fish. 
One  wishes  to  find  a  function  with  modulus  being  a  picture  of  a 
fish,  which  has  a  Fourier  transform  with  modulus  being  a  picture  of 
a  bird.  Or,  in  terms  of  computer  holography,  find  a  phase  function 
to  assign  to  the  image  of  a  fish  so  that  the  hologram  will  look  like 
an  image  of  a  bird.  Figure  4(a)  shows  the  actual  “bird"  and  “fish" 
binary  patterns  used  for  our  experiment.  For  the  first  iteration,  the 
fish  object  was  random  phase  coded,  Fourier  transformed,  and  the 
modulus  of  the  Fourier  transform  was  replaced  with  the  modulus 
of  the  bird  pattern  shown  in  Fig.  4(a).  The  result  was  inverse 
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Fig.  3.  Bird  transforms  into  fish  ("Sky  and  Water”  by  M.  C.  Eschar). 
This  reproduction  was  authorized  by  the  M.  C.  Eschar  Foundation, 
The  Hague,  Holland/G.W.  Breughel. 


(a) 


Fig.  4.  Example  of  spectrum  shaping,  (a)  Bird  hologram  and 
desired  fish  image;  (b)  fish  output  image  alter  random  phase 
coding  of  input;  (c)  output  image  after  seven  Iterations  of  the 
iterative  algorithm. 


Fourier  trsformed,  yielding  the  very  noisy  output  image  shown  in 
Fig.  4(h).  The  iterative  algorithm  was  then  used  for  seven  itera¬ 
tions.  resulting  in  the  improved  image  shown  in  Fig.  4(c).  For  this 
example,  increasing  the  number  of  iterations  resulted  in  a  further 
improvement  of  the  quality  of  the  image;  that  is,  a  Fourier 
transform  pair  was  found  that  more  closely  satisfied  the  constraints 
in  both  domains. 

Spectrum  shaping  is  also  important  in  computer  holography  for 
reducing  quantization  noise.  The  objective  of  computer 
holography16  is  to  synthesize  a  transparency  that  can  modulate  a 
wavefront  according  to  a  calculated  wavefront,  often  correspond¬ 
ing  to  Fourier  coefficients  (or  samples  of  the  Fourier  transform  of 
an  image)  computed  by  the  discrete  Fourier  transform.  Let  F  - 
&  [f]  be  the  desired  wavefront  modulation  and  f  be  the  complex- 
valued  function  describing  the  desired  image.  Due  to  the  limitations 
of  the  recording  devices  and  materials  used  to  synthesize  computer 
holograms,  it  is  often  not  possible  to  represent  exactly  any  arbitrary 
complex  Fourier  coefficient.  An  extreme  example  of  this  is  the 


(b) 


Fig.  5.  Computer-simulated  images  from  kinoform.  (a)  object  ran¬ 
dom  phase  coded;  (b)  alter  eight  iterations  of  the  iterative 
algorithm. 


kinoform,1'  which  allows  nearly  continuous  phase  control  by  vary¬ 
ing  the  thickness  of  the  recording  medium,  but  which  quantizes  the 
modulus  to  a  single  level.  (If  the  gray-level  recording  dev  ice  used  to 
synthesize  a  kinoform  has  a  finite  number  of  gray  levels,  then  the 
phase  is  quantized  as  well.)  The  desired  coefficient  F  is  only  ap¬ 
proximated  by  the  quantized  value  F/  F  .  Since  only  the  squared 
modulus  (the  intensity)  of  the  image  is  observed,  one  is  free  to 
choose  the  phase  of  the  object  (phase  code  the  object)  in  such  a  way 
as  to  reduce  the  variance  (dynamic  range)  of  F  .  In  this  way  the 
quantization  noise  in  kinoforms  and,  to  a  lesser  extent,  in  other 
types  of  computer-generated  holograms  can  be  greatly  reduced. 
Random  phase  and  various  deterministic  phase  codes11'  cause  con¬ 
siderable  reduction  in  the  variance  of  F  ,  but  substantial  errors  re¬ 
main.19 

It  was  for  the  kinoform  application  that  the  iterative  algorithm 
was  first  invented. Figure  5  shows  an  example  of  its  use  for  this 
synthesis  problem.'  Figure  5(a)  shows  the  image  resulting  when  the 
input  image  was  random  phase  coded,  encoded  as  a  kinoform  in 
the  Fourier  plane,  and  reconstructed  by  inverse  Fourier  transfor¬ 
mation.  The  ideal  image  would  be  the  binary  (  0  or  I )  block  letters 

SU.  Figure  5(b)  shows  the  improved  result  after  eight  iterations  of 
the  iterative  algorithm.  In  this  case,  the  image-domain  constraint  is 
that  the  modulus  equal  the  SU  pattern,  and  the  Fourier-domain 
constraint  is  that  the  modulus  equal  a  constant. 

A  problem  very  similar  to  the  kinoform  problem  is  that  ol  syn¬ 
thesizing  a  quasi-random  radar  signal  having  good  autocorrelation 
properties.  Specifically,  one  would  like  to  synthesize  a  radar  signal 
f(t)  which  is  a  pure  phase  function,  i.e.,  f(t)  I ,  over  some  inter¬ 
val  of  time  and  which  has  an  autocorrelation  function  which  ap¬ 
proaches  a  delta-function,  i.e.,  its  Fourier  spectrum  I  ( i  )  :  is  con¬ 
stant  over  the  bandwidth  of  interest.  From  the  examples  shown 
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above,  it  is  obvious  that  the  iterative  method  would  be  an  effective 
tool  for  synthesizing  such  radar  signals. 

Another  spectrum-shaping  application  is  the  phasing  of  elements 
'f  an  array  of  antennas  in  order  to  achieve  a  far-field  pattern  hav¬ 
ing  desirable  properties.  For  example,  one  might  wish  to  phase  the 
antenna  elements  in  such  a  way  as  to  minimize  the  maximum 
sidelobe  of  the  far-field  pattern  or  to  place  nulls  of  the  antenna  pat¬ 
tern  at  several  different  prescribed  locations  simultaneously.  A 
related  application  for  which  the  iterative  method  has  been  used  is 
the  transformation  of  a  Gaussian  laser  beam  into  a  beam  having  a 
more  nearly  rectangular  profile.20 

3.  I.  3.  Ha  ve front  sensing 

The  wavefront  sensing  application  is  very  similar  to  the  electron 
microscopy  problem.  Suppose  that  one  measures  the  image  f(x)  2 
of  a  point  source  using  an  aberrated  optical  system,  where  the  aber¬ 
rations  may  be  due  to  atmospheric  turbulence  or  due  to  the  optical 
system  itself.  Assuming  that  the  aberration  is  a  pure  phase  func¬ 
tion,  then  F(u),  the  Fourier  transform  of  f(x),  has  modulus  F(u) 
equal  to  the  aperture  function  of  the  optical  system.  The  problem  is 
to  reconstruct  the  phase  of  F(u)  given  F(u)  and  f(x)  .  Several  in¬ 
vestigators921-22  have  applied  the  error-reduction  algorithm  to  this 
problem  with  generally  good  results. 

3.2.  Nonnegativity — modulus  constraints 

For  some  reconstruction  problems,  the  physical  quantity  of  interest 
can  be  represented  as  a  nonnegative  function,  and  one  is  able  to 
measure  only  the  modulus  of  its  Fourier  transform  (or  at  least  the 
measured  modulus  information  has  a  much  higher  signal-to-noise 
ratio  than  the  measured  phase).  From  the  Fourier  modulus,  one 
wishes  to  reconstruct  the  Fourier  phase  or,  equivalently,  the  func¬ 
tion  itself.  Since  the  autocorrelation  of  the  function  is  available  as 
the  inverse  Fourier  transform  of  the  squared  Fourier  modulus,21 
this  problem  is  equivalent  to  reconstructing  the  function  from  its 
autocorrelation.  This  problem,  referred  to  as  the  phase  retrieval 
problem  of  optical  coherence  theory,  arises  in  spectroscopy,24  a 
one-dimensional  problem;  in  astronomy,  a  two-dimensional  prob¬ 
lem;  and  in  x-ray  crystallography,2'  a  three-dimensional  problem. 
In  spectroscopy,  the  nonnegative  spectral  density,  g(i<),  is  the 
Fourier  transform  of  the  complex  degree  of  temporal  coherence, 
y(r),  of  which  *>(r)  is  most  easily  measured.  In  x-ray 
crystallography,  the  nonnegative  electron  density  function,  g(x,  y, 
z),  which  is  periodic,  is  the  Fourier  transform  of  the  structure  fac¬ 
tor  Fh|(|,  of  which  Fhli|  is  measured  by  a  diffractometer.  The 
astronomy  problem  will  be  described  in  more  detail  later. 

3. 2. 1.  I  niqueness  of  solutions 

For  the  one-dimensional  problem,  use  of  the  iterative  algorithm  (or 
any  other  method)  to  reconstruct  the  function  from  its  Fourier 
modulus  is  of  limited  interest  since  the  solution  in  the  general  case  is 
usually  not  unique. 2f>-2  The  uniqueness  of  the  solution  for  the  one¬ 
dimensional  problem  can  be  analyzed  using  the  theory  of  analytic- 
functions,  from  which  one  finds  that  additional  solutions  can  be 
generated  by  “flipping  zeros"  of  the  Fourier  transform  analytically 
extended  over  the  complex  plane. 2A-2’  The  additional  "solutions” 
have  the  same  support  as  the  original  function,  but  ate  not 
guaranteed  to  be  nonnegativc;  therefore  one  could  reduce  the 
degree  of  ambiguity  by  generating  all  possible  "solutions"  and  then 
keeping  only  the  nonnegative  ones.2* 

For  certain  special  types  of  one-dimensional  functions,  there  is  a 
high  probability  that  the  solution  is  unique.  For  a  function  having 
two  separated  intervals  of  support,  being  separated  by  an  interval 
over  which  the  function  is  zero,  the  solution  usually  is  unique,29-10 
but  only  if  the  two  intervals  of  support  are  sufficiently  separated. 11 
Another  special  type  of  function  for  which  the  solution  is  usually 
unique  is  one  consisting  of  a  summation  of  a  number  of  delta- 
functions  randomly  distributed  in  space;  for  such  functions,  one 
does  not  need  the  iterative  method —they  can  be  rcconstuctcd  by  a 
simple  noniterative  method  involving  the  product  of  three 


Fig.  6.  Functions  (a)  and  (b)  having  the  same  Fourier  modulus. 


translates  of  the  autocorrelation  function.12 

In  the  event  that  multiple  solutions  do  exist,  it  would  not  appear 
that  the  algorithm  would  be  biased  toward  one  over  another,  and 
one  would  expect  the  algorithm  to  converge  to  different  solutions, 
depending  on  the  initial  input  to  the  algorithm.  For  example.  Fig.  6 
shows  two  functions  having  the  same  Fourier  modulus.  In  a  com¬ 
puter  experiment  using  the  iterative  reconstruction  algorithm  on 
the  functions'  Fourier  modulus,  it  converged  to  one  of  the  solu¬ 
tions  in  about  half  of  the  trials  and  converged  to  the  other  solution 
in  the  other  half  of  the  trials,  depending  on  the  random  number  se¬ 
quences  used  as  the  initial  input  to  the  algorithm. 

For  the  problem  in  two  or  more  dimensions,  it  appears  that  the 
solution  is  usually  unique.  Considering  sampled  functions  defined 
on  a  rectangular  grid  of  points,  Bruck  and  Sodin”  showed  that  the 
existence  of  additional  solutions  is  equivalent  to  the  factorabilitv  of 
a  polynomial  representation  of  the  Fourier  transform.  Since  a 
polynomial  of  one  variable  of  degree  M  can  always  be  factored  into 
M  prime  factors,  there  are  2MI  solutions  in  the  one-dimensional 
case.  Once  again,  only  some  of  the  "solutions"  may  be  non- 
negative.  On  the  other  hand,  polynomials  of  two  or  more  variables 
having  arbitrary  coefficients  arc  only  rarely  factorable;  consequent 
ly,  the  two-dimensional  problem  is  usually  unique.  Attempts  have 
also  been  made  to  extend  this  concept  to  continuous,  as  opposed  to 
discrete,  functions.'4  Although  it  is  always  possible  to  make  up  ex¬ 
amples  in  two  dimensions  that  are  not  unique.”  it  appears  to  be 
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true  that  for  two-dimensional  functions  drawn  from  the  real  world, 
the  solution  is  usually  unique.  The  general  uniqueness  of  the  two- 
dimensional  case  is  indicated  by  experimental  reconstruction  results 
using  the  iterative  algorithm.36  Furthermore,  noise  in  the  Fourier 
modulus  data  has  had  the  effect  of  adding  noise  to  the 
reconstructed  function  rather  than  causing  the  algorithm  to  con¬ 
verge  to  a  radically  different  solution.37 

3.2.2.  Astronomical  reconstruction 

The  problem  of  reconstructing  a  two-dimensional  nonnegative 
function  from  the  modulus  of  its  Fourier  transform  arises  in 
astronomy.  Due  to  atmospheric  turbulence,  the  resolution  at¬ 
tainable  from  large  optical  telescopes  on  earth  is  only  about  one 
second  of  arc,  many  times  worse  than  the  diffraction  limit  imposed 
by  the  diameter  of  the  telescope  aperture.  For  a  five-meter 
telescope  aperture,  the  diffraction-limited  resolution  would  be 
about  0.02  seconds  of  arc — fifty  times  finer.  Despite  atmospheric 
turbulence,  it  is  possible  to  measure  the  modulus  of  the  Fourier 
transform  of  a  space  object  out  to  the  diffraction  limit  of  the 
telescope  using  interferometric  techniques.38'41  The  autocorrelation 
of  the  object  can  be  computed  from  the  Fourier  modulus,  allowing 
the  diameter  of  the  object  to  be  determined.  However,  unless  the 
Fourier  transform  phase  is  also  measured,  it  was  previously  not 
possible  to  determine  the  object  itself,  except  for  some  special 
cases.  Previous  attempts  to  solve  this  problem  had  not  proven  to  be 
practical  for  complicated  two-dimensional  objects. 

The  problem  of  reconstructing  an  object  from  interferometer 
data  can  be  solved  by  the  iterative  method.42,36  The  Fourier- 
domain  constraint  is  that  the  Fourier  modulus  equal  the  Fourier 
modulus  measured  by  an  interferometer,  and  the  function-domain 
constraint  is  that  the  object  function  be  nonnegative.  Figure  7 
shows  an  example.  Fig.  7(a)  shows  a  computer-synthesized  object 
used  for  the  experiment— a  sun-like  disk  having  “solar  flares"  and 
bright  and  dark  “sunspots.”  The  modulus  of  its  Fourier  transform 
is  shown  in  Fig.  7(b).  Figure  7(c)  shows  a  square  of  random 
numbers  used  as  the  initial  input  for  the  iterative  algorithm.  Figures 
7(d),  7(e),  and  7(f)  show  the  reconstruction  results  after  20,  230, 
and  600  iterations,  respectively.  Figure  7(g)  shows  the  initial  input 
for  a  second  trial,  and  the  reconstruction  results  after  2  and  215 
iterations  are  shown  in  Figs.  7(h)  and  7(i),  respectively.  Comparing 
Figs.  7(0  and  7(i)  with  the  original  object  in  Fig.  7(a),  one  sees  that 
for  both  trials,  the  reconstructed  images  match  the  original  object 
very  closely.  Note  that  inverted  solutions  such  as  Fig.  7(0  are  per¬ 
mitted  for  this  problem  since  the  modulus  of  the  Fourier  transform 
of  f(-x)  equals  the  modulus  of  the  Fourier  transform  of  f(x)  for 
real-valued  f(x).  Other  successful  reconstruction  experiments  have 
been  performed  on  data  simulated  to  have  the  types  of  noise  pres¬ 
ent  in  stellar  speckle  interferometry,39  and  it  appears  that  under 
realistic  levels  of  photon  noise  for  fairly  bright  objects,  diffraction- 
limited  images  can  be  reconstructed.37  Initial  expements  have  also 
been  carried  out  on  data  from  telescopes.43 


Fig.  7.  Reconstruction  of  a  nonnegative  (unction  Irom  its  Fourier 
modulus,  (a)  Test  object;  (b)  modulus  ot  its  Fourier  transform;  (c)  ini¬ 
tial  estimate  of  the  object  (first  test);  (d)-(f)  reconstruction  results 
—number  of  iterations:  (d)  20,  (e)  230,  (f)  600;  (g)  initial  estimate  of 
the  object  (second  test);  (h)-(i)  reconstruction  results— number  of 
iterations;  (h)  2,  (i)  215. 


known  finite  extent  (or  support)  and  one  wishes  to  reconstruct  the 
function  with  resolution  appropriate  to  an  aperture  in  the  Fourier 
domain  more  complete  than  the  one  over  which  measurements  were 
actually  taken.  In  some  cases,  the  desired  aperture  is  simply  larger 
than  the  aperture  over  which  measurements  were  taken,  and  so  one 
wishes  to  extrapolate  the  function’s  Fourier  transform,  i.e.,  to  ob¬ 
tain  superresolution  of  the  function.  In  other  cases,  one  has  made 
measurements  over  a  partially  filled  aperture,  in  which  case  one 
wishes  to  interpolate  the  Fourier  transform  of  the  function,  and 
thereby  obtain  an  improved  impulse  response  in  the  function  do¬ 
main. 


3.2.3.  Pupil  reconstruction  and  synthesis 

Another  ease  in  which  one  may  want  to  reconstruct  a  two- 
dimensional  nonnegative  function  from  its  Fourier  modulus  is  in 
pupil  function  determination.  In  a  diffraction-limited  optical 
system,  the  point-spread  function  is  the  squared  Fourier  modulus 
of  the  system’s  pupil  function.  Equivalently,  the  optical  transfer 
function  is  the  autocorrelation  of  the  pupil  function.44  Given  the 
point-spread  function  at  a  given  location  in  an  image  plane,  one 
could  use  the  iterative  algorithm  to  retrieve  the  corresponding  pupil 
function,  in  a  way  that  is  mathematically  equivalent  to  the 
astronomy  problem.  Turning  this  problem  around,  one  could  use 
the  iterative  algorithm  to  synthesize  (design)  a  pupil  function  that 
would  yield  a  given,  desired  point-spread  function  while  possibly 
satisfying  other  desirable  constraints  as  well. 

3.3.  Finite  extent— measurement  over  part  of  an  aperture 

In  a  number  of  reconstruction  problems,  there  is  a  function  of 


3.3./.  Extrapolation  or  superresolution 

The  error-reduction  algorithm  was  first  applied  to  the  extrapolation 
(or  superresolution)  problem  by  Gerchberg.43  Much  has  been  writ¬ 
ten  about  the  iterative  algorithm,  specifically  the  error-reduction 
algorithm,  as  it  relates  to  this  problem,  including  various  ways  of 
understanding  the  algorithm  (see  the  end  of  Sec.  2)  and  proofs  of 
convergence.101213'46  48  For  this  particular  problem,  the  nature  of 
the  constraints  makes  it  possible  to  implement  the  algorithm  by  a 
feedback  optical  processor49,30  taking  on  the  order  of  10  9  seconds 
per  iteration  even  for  the  two-dimensional  ease.  Marks  and  Smith 
describe  these  matters  in  detail  elsewhere  in  this  volume. 

3.3.2.  Interpolation 

In  tomographic  imaging  systems,  many  projections  ot  the  object 
are  measured,  each  projection  yielding  information  about  a  slice 
through  the  Fourier  transform  of  the  object.  W  hen  measurements 
over  only  a  limited  cone  ol  angles  are  made,  the  effective  aperture 
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in  the  Fourier  domain  has  gaps,  and  ihe  impulse  response  of  ihe 
system  is  highly  irregular.  In  applying  the  iterative  algorithm  to  this 
problem,'1'2  the  function-domain  constraint  is  the  finite  extent 
and  nonnegativity  of  the  object,  and  the  Fourier  domain  constraint 
is  that  the  Fourier  transform  equal  the  measured  Fourier  transform 
over  the  measurement  aperture. 

A  problem  similar  to  the  tomography  problem  arises  in  radio 
astronomy.  The  radio  sky  brightness  map  is  a  two-dimensional 
real,  nonnegative  function  which  is  the  Fourier  transform  of  the 
complex  visibility  function.  The  visibility  function  is  measured  by 
radio  interferometry,  and  in  the  case  of  long-baseline  in¬ 
terferometry,  the  visibility  function  is  measured  only  over  a  limited 
set  of  “tracks"  in  the  Fourier  domain,  resulting  in  a  partially-filled 
effective  aperture.  The  error-reduction  algorithm  has  been  used  to 
obtain  improved  maps  by,  in  effect,  interpolating  the  visibility 
function  to  fill  in  the  area  between  the  tracks."  For  this  problem, 
the  constraints  on  the  brightness  map  are  that  it  be  nonnegative  and 
be  zero  outside  the  known  field  of  view.  In  the  visibility  plane,  the 
constraint  is  that  the  complex  visiblity  function  equal  the  measured 
value  within  the  area  of  the  tracks. 

3.4.  Modulus — quantized  values 

As  mentioned  earlier  in  connection  with  spectrum  shaping,  in  com¬ 
puter  holography  one  may  wish  to  encode  the  Fourier  transform  of 
an  image  as  a  computer-generated  hologram,  but  some  types  of 
computer-generated  holograms  can  encode  only  certain  quantized 
complex  v  alues.  The  kinoform  example  discussed  earlier  is  a  special 
type  of  quantization.  A  more  general  example  is  the  Lohmann 
hologram,'4  for  which  the  modulus  and  phase  of  a  complex  sample 
are  determined  by  the  area  and  relative  position,  respectively,  of  an 
aperture  within  a  sampling  cell.  The  number  of  allowable  quantized 
values  is  determined  bv  the  number  of  resolution  elements,  of  the 
recording  device  used  to  fabricate  the  hologram,  used  to  form  one 
cell.  For  this  synthesis  problem,  the  function-domain  constraint  is 
that  the  modulus  of  the  function  equal  the  desired  image  modulus 
and  the  Fourier-domain  constraint  is  that  the  complex  Fourier 
coefficients  fall  on  a  prescribed  set  of  quantized  values.  F.x- 
periments  have  shown  that  synthesizing  such  a  Fourier  transform 
pair  is  possible  using  the  iterative  algorithm.'5-’  For  example.  Fig. 
8(a)  shows  a  simulation  of  an  image  produced  by  a  Lohmann 
hologram  having  only  lour  modulus  and  four  phase  quantization 
levels  when  the  image  was  random  phase  v.  Jed.  Figure  8(b)  shows 
the  image  after  13  iterations,  a  considerable  improvement.  This 
problem  is  one  of  a  more  general  class  of  problems  regarding  the 
transmission  of  coded  data. 

3.5.  Finite  extent— phase 

Finally,  the  iterative  algorithm  has  been  used  to  reconstruct  Ihe 
modulus  of  a  band-limited  signal  from  its  phase. 5h-'  Or,  looking  at 
it  in  another  way,  given  that  a  function  has  finite  extent  and  given 
the  phase  of  its  Fourier  transform,  reconstruct  the  modulus  of  its 
Fourier  transform.  For  this  application,  it  has  been  shown  that  for 
a  wide  class  of  conditions  the  solution  is  unique. 'h  This  application 
will  be  discussed  further  in  Sec.  4. 

4.  ALGORITHM  CONVKRGKNCK  AM) 

A( C KLFRATKD  ALGORITHMS 

V  mentioned  in  Sec.  2,  the  basic  iterative  algorithm  depicted  in 
Fig.  I .  referred  to  as  the  error-reduction  algorithm,  has  been  shown 
to  converge  for  some  applications.  In  this  section,  the  convergence 
is  proven  for  all  applications.  In  addition,  modified  algorithms  that 
often  converge  much  faster  than  the  error-reduction  algorithm  arc 
discussed 

4.1.  (  unvergence  of  the  error-reduction  algorithm 

lor  the  error  reduction  algorithm,  the  mean-squared  error  can  be 
defined  m  general  bv  Fq.  (7)  or  Fq.  (8).  It  is  a  normalized  version 
ot  the  integral  over  the  square  of  the  amount  by  which  the  cotn- 


(b) 


Fig.  8.  Computer-simulated  images  from  hologram  with  four 
magnitude  and  four  phase  quantized  levels,  (a)  Object  random 
phased  coded;  (b)  after  13  iterations  of  the  iterative  method. 


puled  function  (or  the  computed  Fourier  transform)  violate'  the 
constraints  in  the  appropriate  domain.  When  the  mean-squared  er 
ror  is  zero,  then  a  Fourier  transform  pair  has  been  tound  that 
satisfies  all  the  constraints  in  both  domains. 

Consider  again  the  steps  in  the  error-reduction  algorithm 
described  in  Sec.  2.  The  k1"  iteration  starts  with  an  estimate  gk(xt 
that  satisfies  the  function-domain  constraints.  For  any  coordinate, 
x,  the  complex  values  that  gtx)  can  have  that  satisfy  the  function- 
domain  constraints  form  some  set  of  points  in  phasor  space.  I  or 
example,  if  the  modulus  must  equal  fix)  .  then  the  set  ol  such 
points  is  a  circle  of  radius  f(x)  in  phasor  space;  if  the  function 
must  be  nonnegative,  then  the  set  of  such  points  is  the  hall  line  on 
the  nonnegative  real  axis.  The  function  estimate  gk(x)  is  Fourier 
transformed,  yielding  Cik(u).  file  next  step  in  the  algorithm  is  to 
form  Gk(u)  by  changing  Gk(u)  by  the  smallest  possible  amount  that 
allows  it  to  satisfy  the  Fourier-domain  constraints.  Gk(ut  is  then  in¬ 
verse  Fourier  transformed,  yielding  gk(\)  in  the  function  domain 
In  the  final  step,  gk  .  j(x)  is  formed  by  changing  gk(\)  by  the 
smallest  amount  that  allows  it  to  satisfy  the  function-domain  con¬ 
straints.  Now  consider  the  u  inortnalized  squared  error,  given  bv 
the  numerators  in  I  qs.  (7)  and  (8).  In  the  Fourier  domain,  the  tin 
normalized  squared  error  at  the  k'1'  iteration  is 

oc 

cl  k  /  Gk(u)  ( ik(u )  -  du 

X 

(14) 
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j  |gk(x)-8k(x)l2dx, 

-oo 


where  the  second  line  in  this  equation  results  from  Parseval’s 
theorem.  The  unnormalized  squared  error  in  the  function  domain 
at  the  k1*1  iteration  is  given  by 


00 

e0k  =  J  !gk  +  iOO-gk(*)l2dx. 

-oo 


(15) 


Both  g^fx)  and  gk  +  ,(x)  by  definition  satisfy  the  function-domain 
constraints.  Also  at  any  given  coordinate  x,  gk  +  [(x)  is  the  point  in 
phasor  space  satisfying  the  function-domain  constraints  that  is 
closest  to  gk(x).  Therefore,  for  all  values  of  x. 


8k  +  |(*)-8k(*)i  5  i  gk(x)  '  8k(x)  I  •  <16> 

where  equality  holds  only  if  gk(x)  is  just  as  close  in  phasor  space  to 
gk(x)  as  gk  +  |(x)  is.  When  there  is  a  point  in  phasor  space  satisfying 
the  constraints  that  is  closer  to  gk(x)  than  gk(x)  is,  then  the  left- 
hand  side  of  the  expression  above  is  strictly  less  than  the  right-hand 
side.  Therefore,  combining  Eqs.  (14)-(16), 

e0k  s  4k  <17) 

for  a  given  iteration.  From  the  perfect  symmetry  of  the  error- 
reduction  algorithm,  as  seen  from  Fig.  1,  a  similar  result  holds 
when  one  completes  the  iteration  by  satisfying  the  function-domain 
constraints,  thereby  forming  gk  +  |(x),  and  continues  the  next  itera¬ 
tion  by  Fourier  transforming  gk  +  ,(x)  and  causing  its  transform  to 
satisfy  the  Fourier-domain  constraints.  One  then  finds  that 


4.k  +  i  s  «ok  *4k-  (18) 

Therefore,  the  unnormalized  squared  error  can  only  decrease  (or  at 
least  not  increase)  at  each  iteration.  Since  the  normalized  mean- 
squared  error  is  simply  proportional  to  the  unnormalized  squared  er¬ 
ror,  a  similar  result  holds  for  the  errors  defined  by  Eqs.  (7)  and  (8). 

While  the  error-reduction  algorithm  converges  to  a  solution  suf¬ 
ficiently  fast  for  some  applications,  it  is  unbearably  slow  for  others. 
In  most  cases,  the  error  is  reduced  rapidly  for  the  first  few  itera¬ 
tions,  and  then  much  more  slowly  for  later  iterations. 


Fig.  9.  Block  diagram  of  the  system  tor  the  input-output  concept. 


algorithm,  the  input  is  not  necessarily  an  estimate  of  the  function 
or  a  modification  of  the  output,  nor  does  it  have  to  satisfy  the  con¬ 
straints;  instead,  it  is  viewed  as  the  driving  function  for  the  next 
output.  This  viewpoint  allows  one  a  great  deal  of  flexibility  and  in¬ 
ventiveness  in  selecting  the  next  input  and  allows  the  invention  of 
an  algorithm  that  converges  more  rapidly  to  a  solution.  As  will  be 
seen  later,  the  “input-output  algorithm”  actually  comprises  a  few 
different  algorithms,  all  of  which  are  based  on  the  input-output 
point-of-view. 

How  the  input  should  be  changed  in  order  to  drive  the  output  to 
satisfy  the  constraints  depends  on  the  particular  problem  at  hand. 
The  analysis  given  in  the  appendix  for  a  specific  application  can  be 
generalized  as  follows.  Consider  what  happens  when  an  arbitrary 
change  is  made  in  the  input.  Suppose  that  at  the  kth  iteration  the  in¬ 
put  g|((x)  results  in  the  output  gk(x).  Further,  suppose  that  the  input 
is  then  changed  by  adding  Ag(x): 


4.2.  Input-output  algorithms 

Resulting  from  an  investigation  into  the  problem  of  the  slow  con¬ 
vergence  of  the  error-reduction  algorithm,  a  new  and  faster- 
converging  algorithm  was  developed,  the  input-output 
algorithm. ,5-58.7-,M2  xhe  input-output  algorithm  differs  from  the 
error-reduction  algorithm  only  in  the  function-domain  operation. 
The  first  three  operations — Fourier  transforming  g(x),  satisfying 
Fourier  domain  constraints,  and  inverse  Fourier  transforming  the 
result— are  the  same  for  both  algorithms.  Those  three  operations, 
if  grouped  together  as  shown  in  Fig.  9,  can  be  considered  as  a 
nonlinear  system  with  an  input  g(x)  and  an  output  g'(x).  A  prop¬ 
erty  of  this  system  is  that  its  output  is  always  a  function  having  a 
Fourier  transform  that  satisfies  the  Fourier-domain  constraints. 
Therefore,  if  the  output  also  satisfies  the  function-domain  con¬ 
straints,  then  all  the  constraints  are  satisfied  and  it  is  a  solution  to 
the  problem.  It  is  then  necessary  to  determine  how  to  manipulate 
the  input  in  such  a  way  as  to  force  the  output  to  satisfy  the 
ftinct  ion -domain  constraints. 

For  the  error-reduction  algorithm,  the  next  input  g(x)  is  chosen 
to  be  the  current  best  estimate  of  the  function  satisfying  the 
function-domain  constraints.  However,  for  the  input-output 


gk  +  |(x)  =  gk(x)  +  Ag(x) .  (19) 

Then  one  would  expect  the  new  output  resulting  from  gk  4  |(x)  to 
be  of  the  form 

gk  +  ,<x)  =  gk(x)  +  oAg(x)  +  additional  noise.  (20) 


That  is,  the  expected  (or  statistical  mean)  value  of  the  change  of  the 
output,  due  to  the  change  Ag(x)  of  the  input,  is  aAg(x),  a  constant 
times  the  change  of  the  input.  The  system  shown  in  Fig.  9  is  not 
linear;  nevertheless,  small  changes  of  the  input  tend  to  result  in 
similar  changes  of  the  output.  The  expected  value  of  the  change  of 
the  output  can  be  predicted,  but  its  actual  value  cannot  be 
predicted  since  it  has  a  non-zero  variance.  In  the  equation  above, 
this  lack  of  predictability  is  indicated  by  the  “additional  noise" 
term.  The  constant  a  depends  on  the  statistics  of  Gk(u)  and  F(u) 
and  on  the  Fourier-domain  constraints. 

If  the  output  gk(x)  does  not  satisfy  the  function-domain  con¬ 
straints  and  if  gk(x)  +  Agj(x)  does,  then  one  might  tr>  to  drive  the 
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output  to  satisfy  the  constraints  by  changing  the  input  in  such  a 
way  as  to  cause  the  output  to  change  by  Ag^x).  Accordir.  -  the 
equation  above,  the  change  of  the  input  that  will,  on  the  average, 
cause  a  change  Ag^x)  of  the  output  is 

Ag(x)  =  a 'Agjjfx) .  (21) 

Thus  a  logical  choice  for  the  new  input  is 

8k  +  lOO  =  8k(*)  +  tfAg^x) ,  (22) 

where  0  is  a  constant  ideally  equal  to  a'1,  and  where  Agd(x)  is  a 
function  such  that  gk(x)  +  Agjfx)  satisfies  the  function-domain 
constraints.  If  a  is  unknown,  then  a  value  of  0  only  approximately 
equal  to  a'1  will  usually  work  nearly  as  well.  The  use  of  too  small  a 
value  of  0  in  Eq.  (22)  will  only  cause  the  algorithm  to  converge 
more  slowly.  The  noise-like  terms  in  Eq.  (20)  are  kept  to  a 
minimum  by  minimizing  j  /3Agd(x)  | . 

As  mentioned  earlier,  for  the  input-output  algorithm  gk(x)  is  not 
necessarily  an  estimate  of  the  function;  it  is  instead  the  driving 
function  for  the  next  output.  Therefore,  it  does  not  matter  whether 
its  Fourier  transform,  Gk(u),  satisfies  the  Fourier-domain  con¬ 
straints.  Consequently,  for  the  input-output  algorithm,  the  mean- 
squared  error,  E^,  is  unimportant;  Ej  is  the  meaningful  quality 
criterion.  When  computing  Ec  for  the  input-output  algorithm,  the 
gk  +  ,(x)  that  one  should  use  in  the  integrand  of  Eq.  (8)  is  the  one 
determined  by  the  error-reduction  algorithm  rather  than  the  one 
computed  by  the  input-output  algorithm.  That  is,  E0  should  still  be 
a  measure  of  the  amount  by  which  the  output,  gk(x),  violates  the 
constraints. 

Another  interesting  property  of  the  system  shown  in  Fig.  9  is  that 
if  an  output  g'(x)  is  used  as  an  input,  then  its  output  will  be  itself. 
Since  the  Fourier  transform  of  g  (x)  already  satisfies  the  Fourier- 
domain  constraints,  g  (x)  is  unaffected  as  it  goes  through  the 
system.  Therefore,  no  matter  what  input  actually  resulted  in  the 
output  g'(x),  the  output  g  (x)  can  always  be  considered  to  have 
resulted  from  itself  as  an  input  From  this  point  of  view,  another 
logical  choice  for  the  new  input  is 

8k  +  |(x)  =  8k<*)  +  0Agd(x)  (23) 

Note  that  if  0  =  1  in  Eq.  (23),  then  this  version  of  the  input- 
output  algorithm  reduces  to  the  error-reduction  algorithm.  Since 
the  optimum  value  of  0  is  usually  not  unity,  the  error-reduction 
algorithm  can  be  looked  on  as  a  suboptimal  subset  of  one  version 
of  the  more  general  input-output  algorithm.  Depending  on  the 
problem  being  solved,  other  variations  in  Eqs.  (22)  and  (23)  may  be 
successful  ways  for  choosing  the  next  input. 

In  order  to  implement  the  input-output  algorithm  using  Eq.  (22) 
or  (23),  one  chooses  Agd(x)  according  to  the  function-domain  con¬ 
straints.  In  general,  a  logical  choice  is  the  smallest  value  of  Agd(x) 
for  which  gk(x)  +  Agd(x)  satisfies  the  function-domain  constraints. 
At  those  values  of  x  for  which  gk(x)  already  satisfies  the  function- 
domain  constraints,  one  would  set  Agd(x)  =  0.  At  those  values  of  x 
for  which  gk(x)  violates  the  function-domain  constraints,  examples 
of  logical  choices  of  Agd(x)  for  various  applications  are  as  follows. 
For  the  astronomy  problem  and  other  applications  requiring  the 
function  to  be  nonnegative,  choose  Agd(x)  =  -gk(x)  where  gj^(x)  is 
negative.  For  applications  requiring  the  function  to  be  of  finite  ex¬ 
tent,  choose  Agd(x)  =  -gk(x)  for  x  outside  the  known  region  of 
support.  For  applications  requiring  the  function  to  have  modulus 
equal  to  f(x)  ,  choose 

gk(x) 

Agd(x)  =  f(x) - gk(x) .  (24) 

gk(x)' 


In  addition  to  the  values  of  Ag^x)  given  above,  there  are  other 
choices  that  are  successful  when  used  in  Eqs.  (22)  and  (23).  Any 
Agd(x)  that  moves  g(x)  in  the  general  direction  of  satisfying  the 
function-domain  constraints  will  usually  result  in  an  algorithm  that 
works;  suboptimum  choices  of  Agd(x)  and  of  0  in  Eq.  (22)  or  Eq. 
(23)  result  in  algorithms  that  converge  less  rapidly  than  the  op¬ 
timum.  Two  examples  of  other  algorithms  that  converge  more 
rapidly  than  the  “logical”  ones  described  in  the  preceding 
paragraph  are  as  follows.  For  applications  requiring  the  function  to 
have  modulus  equal  to  f(x)  ,  it  was  noticed  that  the  difference  in 
phase  between  gk(x)  and  gk(x)  tends  to  have  the  same  sign  as  the 
change  of  phase  of  g^(x)  from  one  iteration  to  the  next.  In  order  to 
anticipate  the  direction  that  the  phase  's  changing,  one  could 
choose  a  Agd(x)  that  tends  to  rotaif  he  phase  angle  of  the  new 
input  toward  that  of  the  last  output.  That  is,  a  good  choice  for  the 
desired  change  of  the  output  is 


gk(x) 

Agd(x)  =  f(x)  — - -  -  gk(x) 

gk(x) 


gk(x)  gk(x) 

- f(x)  i  — — 

8k(x)  gk(x) 


(25) 


in  which  the  first  component  boosts  (or  shrinks)  the  magnitude  of 
the  output  to  match  f(x)  and  the  second  component  rotates  the 
phase  angle  of  the  input  toward  the  phase  angle  of  the  output.  For 
the  astronomy  problem,  it  was  found  that  a  particularly  successful 
algorithm  was  to  use  Eq.  (23)  at  those  points  where  the  constraints 
were  satisfied  and  use  Eq.  (22)  at  those  points  where  the  constraints 
were  violated,  i.e. , 


|gk(x),  where  constraints  satisfied 

(26) 

gk(x)  -  0gk(x),  where  constraints  violated 

Furthermore,  it  was  found  that  even  faster  convergence  can  be  ob¬ 
tained  by  alternating  between  the  above  equation  and  the  error- 
reduction  algorithm  every  few  iterations. 

Unlike  the  error-reduction  algorithm,  the  input-output  algorithm 
is  not  guaranteed  to  converge;  in  fact  the  error  may  even  increase 
for  some  of  the  iterations.  However,  the  input-output  algorithm  is 
much  less  prone  to  stagnation  and  therefore  in  practice  converges 
much  faster  than  the  error-reduction  algorithm.  In  some  instances 
during  the  input-output  iterations,  Ec  may  even  increase  although 
the  visual  appearance  of  the  image  improves.  This  behavior,  which 
is  poorly  understood,  is  described  further  in  Ref.  59. 

From  the  paragraphs  above,  it  is  seen  that  the  “input-output 
algorithm”  is  really  a  family  of  algorithms.  The  input-output  ap¬ 
proach  is  one  that  can  lead  to  a  number  of  different  algorithms 
based  on  the  manner  in  which  the  nonlinear  system  of  Fig.  9 
behaves.  One  would  hope  that  the  principles  of  control  theory  and 
possibly  other  disciplines  could  be  used  to  shed  further  light  on  this 
system  and  help  to  arrive  at  algorithms  with  still  more  rapid  con¬ 
vergence. 

It  should  also  be  noted  that,  unlike  the  error-reduction 
algorithm,  the  input-output  algorithm  does  not  treat  the  two  do¬ 
mains  in  a  symmetric  manner.  By  reversing  the  roles  of  the  two 
domains,  one  can  arrive  at  a  different  and  possibly  more  advan¬ 
tageous  algorithm. 

4.3.  Relaxation-parameter  algorithm 

A  second  method  of  improved  convergence  is  the  use  of  a  relaxa¬ 
tion  parameter.  In  solving  the  problem  of  reconstructing  the 
magnitude  of  a  band-limited  function  from  its  phase  (or, 
equivalently,  reconstructing  a  function  of  finite  extent  from  the 
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Fig.  10.  Block  diagram  of  the  error-reduction  algorithm  modified  to 
include  a  relaxation  step. 


phase  of  its  Fourier  transform),  Oppenheim,  Hayes,  and  Lim57 
modified  the  error-reduction  algorithm  (Fig.  1)  by  adding  a  relaxa¬ 
tion  step,  as  shown  in  Fig.  10.  Here  the  band-limited  function  is 
taken  to  be  in  the  Fourier  domain.  The  function  g(x)  then  must  be 
of  finite  extent  according  to  the  bandwidth  of  the  Fourier-domain 
function.  In  the  relaxation  step,  gk(x)  is  formed  from  gk(x)  accord¬ 
ing  to 

g£(x)  =  (1  -  tjk)gk.,(x)  +  »ikgk(x) ,  (27) 

and  then  the  new  estimate  gk  +  ((x)  is  formed  from  gk(x)  by  making 
it  satisfy  the  function-domain  constraints.  The  parameter  ijk,  which 
is  a  constant  that  may  vary  from  one  iteration  to  the  next,  is  the 
relaxation  parameter.  For  t;k  =  l,gk(x)  =  gk(x)  and  this  reduces  to 
the  error-reduction  approach.  For  rjk  =  0,  gk(x)  =  gk_,(x),  that  is, 
the  result  from  the  previous  iteration  is  used.  Other  values  of  ijk 
give  a  linear  combination  of  gk'_,(x)  and  gk(x).  For  the  reconstruc¬ 
tion  of  a  function  of  finite  extent  from  the  phase  of  its  Fourier 
transform  or  from  a  segment  of  its  Fourier  transform  (i.e.,  the 
superresolution  problem),  if  g,'(x)  and  g'(x)  both  satisfy  the 
Fourier-domain  constraint,  then  the  linear  combination  qg[(x)  + 
(1  -  Tj)gj(x)  also  satisfies  the  constraint  in  the  Fourier  domain,  it 
follows  from  this  that  gk(x)  given  by  Eq.  (27)  also  satisfies  the 
Fourier-domain  constraint.  In  those  cases,  it  can  be  shown  that  the 
algorithm  converges  for  0  <  r?k  <  I .  However,  for  other  sets  of 
constraints,  for  example,  given  (*•£  modulus  of  the  Fourier  trans¬ 
form,  gk(x)  given  by  the  equation  above  does  not  generally  satisfy 
the  Fourier-domain  constraints  and  so  the  relaxation  method  does 
not  strictly  apply. 

The  optimum  value  of  i7k  can  be  determined  as  follows.  Define 
the  function-domain  squared  error  after  the  relaxation  step  as 


The  computation  of  the  relaxation  parameter  by  Eq.  (29)  lakes 
much  less  time  than  the  computation  of  one  (fast)  Fourier 
transform,  and  so  it  does  not  significantly  increase  the  total  com¬ 
putation  time  of  a  single  iteration. 

Use  of  the  relaxation  step  for  the  problem  of  reconstructing  a 
band-limited  function  from  its  phase  resulted  in  an  order  of 
magnitude  improvement  in  the  speed  of  convergence  of  the 
algorithm  over  that  of  the  error-reduction  algorithm/7 

The  relaxation  step  described  above  incorporates  the  optimum 
combination  of  the  current  output  with  the  previous  output.  It  is 
also  possible  to  extend  this  concept  to  include  a  number  of  previous 
outputs,57  which  may  result  in  still  more  rapid  convergence. 

It  should  be  noted  that  the  majority  of  the  work  referenced  in 
Sec.  3  made  use  of  only  the  error-reduction  algorithm.  Improved 
speed  of  convergence  could  be  expected  if  one  of  the  two  ac¬ 
celerated  algorithms  discussed  above  were  employed. 

S.  SUMMARY  AND  COMMENTS 

The  iterative  error-reduction  algorithm,  an  extension  of  the 
Gerchberg-Saxton  algorithm  to  include  various  types  of  con¬ 
straints,  has  been  found  to  be  capable  of  solving  a  wide  range  of 
difficult  problems  in  optics  and  other  fields.  It  can  be  applied  to  the 
reconstruction  of  a  function  (an  object,  wavefront,  signal,  etc.) 
when  only  partial  information  is  available  in  each  of  two  domains, 
or  to  the  synthesis  of  a  function  (wavefront,  signal,  etc.)  having 
desired  properties  in  each  of  two  domains.  The  iterative  algorithm 
is  reasonably  fast  for  most  applications,  since  the  major  computa¬ 
tional  burden,  two  Fourier  transforms  per  iteration,  can  be  ac¬ 
complished  using  the  fast  Fourier  transform  (FFT)  algorithm.  The 
iterative  algorithm  has  been  shown  to  outperform  alternative 
methods  of  solving  these  classes  of  problems  both  because  of  its 
speed  and  its  tolerance  of  noise.4'9  For  some  applications,  a  large 
number  of  iterations  is  required  for  convergence  of  the  error- 
reduction  algorithm.  This  situation  can  be  remedied  by  using  an 
algorithm  with  accelerated  convergence,  such  as  the  input-output 
algorithm  or  an  algorithm  employing  a  relaxation  step. 

The  iterative  algorithm  has  been  in  use  for  only  a  few  years,  yet  it 
has  already  found  numerous  applications;  and  methods  of  improv¬ 
ing  the  algorithm  have  been  devised.  Nevertheless,  it  is  safe  to 
predict  that  it  will  be  used  in  the  future  to  solve  new  problems  not 
discussed  here,  and  it  is  hoped  that  further  improvements  of  the 
algorithm  will  be  discovered. 

POSTSCRIPT 

As  this  book  goes  to  print,  further  developments  relating  to  the 
iterative  algorithm  are  occurring  at  a  rapid  pace.  It  has  been  un¬ 
covered  that  an  algorithm  equivalent  to  Gerchberg's4'  error- 
reduction  algorithm  for  extrapolation  was  proposed  by  Ville60  in 
1956,  although  approached  from  a  different  point  of  view.  Rela¬ 
tionships  between  ihe  error-reduction  algorithm  and  gradient 
search  methods  have  been  discovered59'61  '6:  and  uncovered.6’  And 
further  work  on  various  applications  is  being  reported.64  *' 


/  gk(x)  2  dx  , 
7 


(28) 


where  the  region  of  integration,  >,  is  the  region  over  which  the 
function  is  known  to  be  zero.  Setting  equal  to  zero  the  derivative  of 
e(-,  with  respect  to  >jk,  and  solving  for  ijk,  one  finds  the  optimum 
value  of  ijk  to  be  given  by 


7k 


-Re 


/ 


gk.,(x)|gk(x)  -  gk.,(x)|*  dx 


/ 


gk(x) 


gk.|(x)  ’  dx 


(29) 


APPENDIX:  ANALYSIS  OF  THE  INPUT-OUTPUT 
SYSTEM 

Consider  the  synthesis  problem  for  kinoforms,  for  which  the 
Fourier  modulus  is  set  equal  to  a  constant.  Suppose  that  the  input 
g(x)  to  a  kinoform  system  results  in  the  output  g  (x).  The  kinoform 
has  a  transmittance  G  (u)  =  K  cxp[i<Mu)|,  where  O(u)  is  the  phase 
of  G(u)  =  G(u)  exp(i<s(u))  =  [g(x)j,  and  K  is  a  constant.  The 

resulting  image  is  g  (x)  '((;  (u)].  Now  consider  what  happens 

when  a  change  Ag(x)  is  made  in  the  input.  As  illustrated  in  the 
phasor  diagrams  in  Fig.  Al ,  the  change  Ag(x)  of  the  input  causes  a 
change  AG(u)  of  its  Fourier  transform,  which  causes  a  change 
AG  (u)  of  the  kinoform  and  a  corresponding  change  Ag  ( x) 
'|AG  (u)|  of  the  output  image.  The  goal  here  is  to  determine 
the  relationship  between  the  change  Ag  (x)  of  the  output  and  the 
change  Ag(\)  of  the  input.  Figure  A2  shows  the  relationship  bc- 


SPIE  Vo /  373  Transformations  in  Optica!  Signal  Processing  ft  981 1  157 


Fig.  A1.  A  change  Ag  of  the  input  results  in  a  change  AG'  of  the 
klnoform  and  a  change  of  Ag '  of  the  output. 


Fig.  A2.  Relationship  between  aG',  the  change  of  the  klnoform, 
and  two  components  of  AG,  the  Fourier  transform  of  the  change  of 
the  Input. 


tween  AG'(u)  and  two  orthogonal  components  of  AG(u).  By 
similar  triangles,  for  AG  <<  G  , 

AG '  (u)  =  AG'(ii) - - - ,  ( A 1 ) 

G(u) 

where  the  two  orthogonal  components  of  AG(u)  are 

AGr(u)  =  AG(u)  cos /3(u)  e‘*<u'  (A2) 

parallel  to  G(u),  and 

AG'(u)  =  AG(u)  sintf(u)eil0,u>tlr/2l  (A3) 

orthogonal  to  G(u);  and 

AG(u)  =  AGr(u)  +  AG'(u)  =  AG(u)  eil<6(u)  +  ,j<u)l,  (A4) 

where  d(u)  is  the  angle  between  AG(u)  and  G(u).  Only  one  of  the 
two  orthogonal  components  of  AG(u),  namely  AG’(u),  contributes 
to  AG’(u). 
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In  order  to  compute  the  expected  change  of  the  output, 
EfAg'(x)],  treat  the  phase  angles  /3(u)  and  the  magnitudes  |  G(u)  |  as 
random  variables.  Inserting  |AG(u)|  from  Eq.  (A4)  into  Eq.  (A3), 
one  obtains 

AG'(u)  =  AG(u)  e-'l^u) +  ^<u)l  sin  0(u)  ei<s<u>  e'*/2 


=  AG(u)[sin2/3(u)  +  i  sin  0(u)  cos  y3(u)]  . 


For  0(u)  uniformly  distributed  over  [0,  2r],19  the  expected  value  of 
AGl(u)  is 

E[AG‘(u)1  =  AG(u)  +  i-oj  =  AG(u) .  (A6) 

Therefore,  the  expected  value  of  the  change  of  the  output  is,  using 
Eqs.  (Al)  and  (A6)  and  assuming  that  the  magnitudes  G(u)j  are 
identically  distributed  random  variables19  independent  of  <3(u), 


E[Ag  (x)J  =  E  ^  9 (AG') 


9  [E(AG'))  =  9  E(AG')  E 


=  9  —  AG(u)  E  (—  \  =  —  Ag(x)E[  — 

.  2  J  \  iC  /  2  y  G 


That  is,  the  expected  change  of  the  output  is  n  times  the  change  of 
the  input,  giving  us  the  second  term  in  Eq.  (20),  where  a  = 
(l/2)E(K/  G1).  After  a  few  iterations,  G(u)  will  not  differ 
greatly  from  K;  then  a  -  1/2. 

Similarly,  the  variance  of  the  change  of  the  output  can  be  shown 
to  be5* 

E(  Ag'(x)  :]  -  E[Ag  (x)]  2 


Ag(x  ' )  -  dx  . 

(A8) 


where  A  is  the  area  of  the  image.  That  is,  the  variance  of  the  change 
of  the  output  Ag'(x)  at  any  given  x  is  proportional  to  the  integrated 
squared  change  of  the  entire  input.  The  predictability  of  Ag  (x), 
and  the  degree  of  control  with  which  one  can  manipulate  it, 
decreases  as  one  makes  larger  changes  of  the  input.  The  difference 
between  the  actual  change  of  the  output  and  the  expected  change  of 
the  output  given  by  Eq.  (A7)  is  what  is  meant  by  the  additional 
noise  term  in  Eq.  (20).  If,  after  a  few  iterations,  G(u)  -  k.  then 
in  Eq.  (A8)  the  factor  (l/4)|2E(K2/  G  2)  -  [E(K  G  )]2I  =  14. 

Equations  (A7)  and  (A8)  are  a  justification  for  the  input-output 
concept:  small  changes  of  the  input  result  in  similar  changes  of  the 
output,  and  so  the  output  can  be  driven  to  satisfy  the  constraints  b> 
appropriate  changes  of  the  input,  as  in  Eqs.  (22)  and  (23). 
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A  simple  recursive  algorithm  is  proposed  for  reconstructing  certain  classes  of  two-dimensional  objects  from  their 
autocorrelation  functions  (or  equivalently  from  the  modulus  of  their  Fourier  transforms — the  phase-retrieval  prob¬ 
lem).  The  solution  is  shown  to  be  unique  in  some  cases.  The  objects  contain  reference  points  not  satisfying  the 
holography  condition  but  satisfying  weaker  conditions.  Included  are  objects  described  by  Fiddy  et  at.  |Opt.  Lett. 
8, 96  (1983)|  satisfying  Eisenstein's  thorem. 


INTRODUCTION 

In  a  number  of  disciplines,  including  astronomy,  x-ray  crys¬ 
tallography,  electron  microscopy,  and  wave-front  sensing,  one 
encounters  the  phase-retrieval  problem.  One  wishes  to 
reconstruct  f(m,  n),  an  object  function,  from  |F(p,  q)|,  the 
modulus  of  its  Fourier  transform,  where 

F{p,q)  =  |F(p,  t?)|exp[t^(p,  <?)]  =  n)] 

At—  1  N-l 

=  L  L  /(m,  n)exp(-i2ir(mp/M  +  nq/N)],  (1) 

m  =0  rt  *0 

where  m,p  =  0, 1 . M  -  1  and  n,q  =  0, 1,  . . .  ,N  -  1.  The 

discrete  transform  is  employed  here  since  in  practice  one  deals 
with  sampled  data  in  a  computer.  The  problem  of  recon¬ 
structing  the  object  from  its  Fourier  modulus  is  equivalent  to 
reconstructing  the  Fourier  phase,  \p(p,  q),  from  the  Fourier 
modulus;  since  once  one  has  the  phase  as  well  as  the  modulus, 
one  can  easily  compute  f(m,  n)  by  the  inverse  (discrete) 
Fourier  transform.  ry(m,  n),  the  (aperiodic)  autocorrelation 
of  f(m.  n),  is  given  by1 

M - 1  N-l 

r/(m,  n)  =  £  £  f(j,k)f*(j  -  m,k  -  n)  (2) 

j  =o  k=  o 

=  J-'\\F(p,q)\z],  (3) 

where  the  asterisk  denotes  complex  conjugate.  Note  that  the 
autocorrelation  is  Hermitian;  ry(— m,  —  n)  =  r/*(m,  n).  Note 
also  that  in  order  to  avoid  aliasing  during  the  computation  of 
|F(p,  q)|-,  it  is  necessary  to  have/(m,  n)  =  0  for  M/2  <  m  < 
M  -  1  and  for  N/2  <  n  <  N  —  1;  this  will  be  assumed 
throughout  this  paper.  Then  there  is  no  difference  between 
the  periodic  (cyclic)  and  aperiodic  autocorrelation.  (For  x-ray 
crystallography  this  is  usually  not  the  case,  and  the  results  of 
this  paper  do  not  apply.)  Since  the  autocorrelation  function 
is  easily  computed  from  the  Fourier  modulus  by  Eq.  (3),  the 
phase-retrieval  problem  is  equivalent  to  reconstructing  an 
object  from  its  autocorrelation  function. 

Several  phase-retrieval  algorithms  have  been  proposed,  all 
of  them  requiring  some  additional  measurements  or  con¬ 
straints  on  the  solution.  Examples  include  a  reference  point 
at  least  one  object  diameter  from  the  object2  (giving  rise  to  the 
holography  condition11),  a  second  intensity  measurement  in 
another  plane4  "r’  (in  electron  microscopy  or  wave-front  sens¬ 


ing),  nonnegativity  and  limited  spatial  extent15  s (in  astrono¬ 
my),  atomic  models9  (in  x-ray  crystallography),  and  objects 
consisting  of  collections  of  points  having  nonredundant 
spacings.10 

Here  it  is  pertinent  to  review  the  case  of  holography. 
Suppose  that  f(m,n)  consists  of  an  object  of  interest,  g(m ,  n ), 
plus  an  unresolved  (delta-function-like)  point,  referred  to  as 
the  reference  point,  i.e., 

f(m,  n)  =  Ab(m  -  m0,  n  -  n0)  +  g(m.  n),  (4) 

where  S(m,  n)  is  a  two-dimensional  (2-D)  Kronecker  delta 
function.  Then  the  autocorrelation  can  be  written  as  the  sum 
of  four  terms, 

ry(m,  n)  =  |A|  25(m,  n)  +  rK(m,n)  +  Ag*(m0  —  m,  tin  -  n ) 
+  A*g(m  +  m0,  n  +  no),  (5) 

the  final  term  of  which  is  the  cross-correlation  of  the  reference 
point  with  the  object  of  interest  and  is  simply  proportional  to 
a  translate  of  the  object  of  interest.  If  the  distance  from  the 
reference  point  to  the  object  of  interest  exceeds  the  diameter 
of  the  object  of  interest,  then  the  fourth  term  in  Eq.  (5)  is 
nonoverlapping  with  the  other  terms,  and  the  object  of  interest 
is  reconstructed  by  simple  inspection  of  the  autocorrelation. 
Then  the  holography  condition  is  satisfied.2-4  If  the  ampli¬ 
tude  and  position  of  the  reference  point  are  unknown  (except 
that  the  reference  point  satisfies  the  holography  condition), 
then  the  object  can  be  reconstructed  only  to  within  a  complex 
factor  A*  and  to  within  a  translation,  and  there  would  be  a 
twofold  ambiguity  as  to  whether  the  object  is  given  by  the 
fourth  term  or  the  third  term  (the  conjugate  image)  of  Eq. 
(5). 

In  this  paper  we  describe  an  algorithm  for  reconstructing 
certain  objects  having  reference  points  that  do  not  satisfy  the 
holography  condition.  For  these  cases  the  reference  points 
may  be  referred  to  as  latent  reference  points,  because  they  do 
not  immediately  yield  the  object  as  would  a  holographic  ref¬ 
erence  point;  rather,  a  degree  of  development  is  required  be¬ 
fore  their  usefulness  emerges. 

In  Section  2  the  question  of  the  uniqueness  of  the  solution 
is  reviewed.  In  Section  3  the  new  reconstruction  algorithm 
is  described  as  it  is  applied  to  three  different  classes  of  objects. 
Additional  comments  on  the  reconstruction  algorithm  are 
included  in  Section  4. 
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2.  UNIQUENESS  OF  THE  SOLUTION 

When  one  measures  only  the  Fourier  modulus,  then  the 
uniqueness  of  the  solution  is  a  central  question.  One  of  course 
always  has  the  twofold  (180°  rotated  or  conjugate  image) 
ambiguity  since  |  ,7\f(m.  n  )||  =  |./'[/*(-m,  -n  )]|; and  trans¬ 
lations  of  f(m,n)  and  the  multiplication  of  f(m.rt)  by  a  con¬ 
stant  phase  factor  exp(tfl)  (where  0  is  a  real  constant)  also  have 
no  effect  on  |F(p,(/l|  If  these  are  the  only  ambiguities,  then 
we  consider  the  solution  of  the  phase-retrieval  problem  to  be 
unique. 

Bruck  and  Sodin11  considered  objects  consisting  of  a  rec¬ 
tangular  grid  of  delta  functions  having  various  complex  am¬ 
plitudes  (or  equivalently,  a  2-D  sequence),  which  have  Fourier 
transforms  that  can  be  expressed  as  polynomials.  These  are 
the  types  of  objects  assumed  by  Eqs.  (1 )  and  (2),  and  we  refer 
to  such  objects  as  sampled  objects.  They  showed  that,  for 
sampled  objects,  a  lack  of  uniqueness  of  the  solution  to  the 
phase-retrieval  problem  is  equivalent  to  the  factorability  of 
the  polynomial,  and  therefore  one-dimensional  (1-D)  objects 
of  length  L  have  a  2/  _) -fold  ambiguity."  This  result  corre¬ 
sponds  to  the  analogous  theory  for  1-D  continuous  functions.12 
On  the  other  hand,  polynomials  of  two  (or  more)  variables  are 
known  to  be  only  rarely  factorable  (i.e.,  they  are  usually  irre¬ 
ducible).  Consequently,  for  2-D  sampled  objects  the  solution 
to  the  phase-retrieval  problem  is  usually  unique.  An  analo¬ 
gous  theory  for  2-D  continuous  functions  is  not  yet  avail¬ 
able. 

Uniqueness  Condition  Due  to  Eisenstein’s  Theorem 

Although  most  2-D  sampled  objects  are,  as  discussed  above, 
uniquely  related  to  the  modulus  of  their  Fourier  transforms, 
it  is  of  interest  to  know  conditions  that  ensure  uniqueness. 
Such  a  condition  was  recently  put  forward  by  Fiddy  et  al. 13 
They  considered  the  class  of  sampled  objects  whose  support 
is  contained  in  the  union  of  a  rectangle  and  an  isolated  point 
(A)  below  and  to  the  right  of  the  rectangle,  as  shown  in  Fig. 
1(a).  By  way  of  example,  the  rectangular  region  in  Fig.  1(a) 
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Kig.  I  Kiddy  Brume-.  Itaintv1  ’  nliiect  la)  FBI)  object  support 
having  two  reference  point*,  .4  and  H.  i  hi  object  support  assumed;  (c> 
autocorrelation  support  The  obiect  is  uniquely  reconstructed  from 
it-  autocorrelation  (unction. 


contains  five  columns  and  four  rows  of  points.  The  object 
must  also  be  nonzero  both  at  point  A  and  at  point  B  in  the 
lower  left  corner  of  the  rectangle.  Points  A  and  B  are  referred 
to  as  the  reference  points,  and  they  do  not  satisfy  the  holog¬ 
raphy  condition.  If  these  conditions  are  satisfied,  then  the 
Fourier  transform  of  the  object  satisfies  Eisenstein’s  theorem, 
making  it  an  irreducible  2-D  polynomial  and  guaranteeing 
that  the  solution  to  the  phase  retrieval  problem  is  unique 
They  demonstrated  the  power  of  these  conditions  by  recon¬ 
struction  experiments  using  the  input-output  iterative  Fou 
rier-transform  algorithm.*1-7  First,  they  performed  a  recon¬ 
struction  experiment  on  the  Fourier  modulus  of  a  particular 
object  that  did  not  have  a  reference  point  A.  After  2 oil  iter¬ 
ations,  a  poor  reconstruction  resulted.  But  when  a  new  object 
was  formed  by  adding  a  reference  point  A  off  its  corner  making 
it  satisfy  the  conditions,  then  a  good  reconstruction  was  oh 
tained  after  only  20  iterations.1 1  Note  that  this  does  not  prove 
that  the  original  object  (without  the  point  A I  was  nonunique: 
the  failure  of  the  iterative  reconstruction  algorithm  may  only 
be  an  indication  of  local  minima  in  the  error  function.  In  fact, 
when  the  reference  point  A  had  a  small  value,  a  poor  recon¬ 
struction  was  obtained  in  spite  of  the  fact  that  irreducibilin 
(and  uniqueness)  was  ensured.  Only  when  a  large  value  lor 
A  was  used  did  the  reconstruction  become  easier  1  Appar¬ 
ently  the  use  of  a  large  enough  value  for  A  also  ensures  that 
there  are  no  local  minima. 

3.  NEW  RECONSTRUCTION  ALGORITHM 

For  certain  classes  of  sampled  objects  having  reference  points 
not  satisfying  the  holography  condition,  we  present  a  new 
reconstruction  algorithm  having  a  fixed  number  of  step* 
This  new  algorithm  is  related  to  the  Dallas''  recursive  algo¬ 
rithm  for  phase  retrieval  from  two  intensity  measurements 
but  requiring  only  a  single  intensity  measurement  (the  Fourier 
modulus)  and  solving  the  equations  in  a  certain  order  such 
that  the  problem  of  a  growing  tree  of  solutions  ’  is  avoided. 
First  the  algorithm  will  be  described  for  the  type  of  object 
described  above,  and  later  for  a  wider  class  of  objects. 

A.  Fiddy-Brames-Dainty  Objects 

For  mathematical  simplicity,  consider  a  sampled  object  whose 
support  is  contained  in  the  regions  shown  in  Fig.  lib).  Its 
uniqueness  properties  are  the  same  as  those  of  the  objects 
considered  in  Fig.  1(a)  since  the  supports  are  mirror  images 
of  one  another.  The  object  can  be  expressed  as  in  Eq.  1 4 1  wit  h 

m(i  =  rtn  =  0: 

fim,  n)  =  Afiim,  n  I  +  g(m.  n  I, 

where g(m.  rt )  is  that  part  of /(m,  n  I  contained  in  the  rectan¬ 
gular  region  of  support,  and  A  =  /(d,  01  *  0.  In  this  case.  elm. 
,i )  is  zero  outside  1  <  m  <  ■)  and  1  <  n  <  K  :  and  it  is  assumed 
that  f(-J.  1 )  =  fj[J,  1 )  =  R  ^  0,  and  glm.  K I  ^  0  for  at  least  one 
value  of  m.  We  will  refer  to  objects  satisfying  these  con¬ 
straints  as  Fidd  v-Brames-Dainty  I  FBI))  objects  having  FBI) 
regions  of  support. 

The  autocorrelation,  r/(m,  n ).  of  fim,  rt )  is  given  by  the  four 
terms  of  Eq.  (5)  with  m„  =  n,,  =  0,  the  supports  of  which  are 
contained  in  the  sets  of  points  illustrated  in  Fig.  lie).  From 
this  figure,  it  can  be  clearly  seen  that  the  rightmost  column 
and  the  uppermost  row  of  r((m ,  n )  are  simply  equal  to  A *e <  m , 
n  V. 
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r/(J,  n )  =  A*g{J,  n), 


(6) 


r/(m,  K)  =  A*g(m,  K),  m  =  l,...,J.  (7) 

Therefore,  for  m  =  J  and  for  n  =  K,  one  can  reconstruct g(m , 
n )  to  within  a  constant  factor  A*  by  simple  inspection  of  ry(m , 
n )  In  effect,  the  holography  condition  is  in  force  for  the  row 
and  column  opposite  reference  point  A,  and  that  row  and  that 
column  are  reconstructed  bv  using  reference  point  A. 

The  value  of  A  can  be  obtained  as  follows:  From  Eq.  (2), 
it  is  seen  that  there  is  only  one  nonzero  term  in  the  summation 
for  the  upper  left  corner  point  in  the  autocorrelation: 


r,(-J  +  l,K  -  1)  =  g(l.K)g*(J,  1)  =  B*gtl,K).  (8) 


Also,  from  Eqs.  (6)  and  (7), 


rf(J,l)  =  A*g(J,\)  =  A*B,  (9) 

ry(l.K)  =  4*g(l,K).  (10) 

Combining  Eqs.  (8)— ( 10)  yields  'ssuming  that  r/t-J  +  1,  K 

-1)^0, 


Ml'  = 


r,[J.  1  )r,*(  1.  K) 


•(-■J  +  1.  K  -  1) 


(11) 


Since  without  loss  of  generality  we  can  arbitrarily  fix  the  phase 
of  any  one  point  in  f(m,n ),  we  set  the  phase  of  A  equal  to  zero; 
A  is  then  gi".  t .  unambiguously  by  the  positive  square  root  of 
Eq.  (II).  II  ■, \~J  +  1 .  K  -  11=0.  then  one  can  obtain  a 
similar  expression  for  |A|-  using  the  first  nonzero  point,  ry(m, 
K  -  II.  to  the  right  of  ry< -J  +  1 ,  K  -  1 >.  Since  A  is  known, 
g(-J.  n  I  and  g(m.  K)  can  be  determined  unambiguously  from 
Kqs.  (6)  and  (7).  Note  that  B  =  gtJ,  1 )  =  ry|«7,  1  )/A  *. 

Having  the  values  of  the  top  row  and  rightmost  column  of 
gtm.  n),  one  can  then  solve  for  the  leftmost  column  in  the 
second  step  of  the  algorithm.  From  Eq.  (2),  the  point  of  the 
autocorrelation  just  below  r/(— J  +  1,  K  -  1)  has  only  two 
nonzero  terms. 


r,<— 7  +  1,  K  -  2)  =  g(l,K)g*(J,  2)  +  g(l,  K  -  1  )g*(J.  1). 

(121 


Solving. 

gtl.K  -  1)  =  |r/<— 7  +  1,  K  -  2)  -  g(l,  K)g*(J ,  2)\/B*. 

<i:d 

where  gt-J,  1 1  =  B  Since  all  the  quantities  of  the  right-hand 
side  of  Eq.  ( 18)  are  known  and  B  ^  0.  one  can  unambiguously 
compute  g 1 1,  K  -  1).  Similarly,  the  next  lower  point  in  the 
autocorrelation  is  given  by 

r/t—J  +  1 ,  K  -  :))  -  gt\.K)g*(-J.  :i)  +  g(l  ,K  -  \)g*tJ.2) 

+  gtl.fi  -  2)g*tJ.  1 1  (141 

Since  all  the  quantities  in  this  linear  equation  are  known  ex¬ 
cept  for  g(  1,  K  -  2),  and  since gtJ,  1 )  ^  0,  one  can  solve  un¬ 
ambiguously  for  gtl.  K  -  2).  In  a  similar  fashion,  one  can 
recursively  solve  for  all  the  values  gt  1,  n )  (the  first  column  on 
the  left)  using  the  values  of  r/(  —  J  +  1.  n  —  1)  in  this  second 
step  of  the  reconstruction.  In  a  sense  the  column  m  =  1  was 
solved  using  the  latent  reference  point  B,  which  required  the 
solution  of  column  m  =  J  before  it  could  become  effective. 

Having  the  first  column  on  the  left  and  the  first  column  on 
the  right  of#(m,  n ),  one  can  then  solve  for  the  second  column 
on  the  right  in  the  third  step,  using  A  as  the  latent  reference 


point.  From  Eq.  (2),  the  points  of  the  autocorrelation  in 
column  (J  -  1)  are  given  by 

,  K 

r,(J  -  \,n)=g(J- l.n)A*+  £  g(J .  k)g*(l .  k  -  n). 

k-n+  l 

(15) 

for  n  =  1 . K  -  1.  Since,  for  any  n,g{J  -  1,  n )  is  the  only 

unknown  in  Eq.  (15),  and  since  A  ^  O.glJ  -  1,  n)  is  uniquely 
determined  from  Eq.  (15).  Thus  the  values  of  g(m.  n)  in 
column  (J  -  1 )  are  reconstructed  using  the  values  in  column 
(J  -  1)  of  the  autocorrelation. 

The  reconstruction  algorithm  continues  in  the  manner 
described  above.  In  the  fourth  step,  one  can  recursively  solve 
for  g(2,  n)  using  the  latent  reference  point  B  and  the  values 
of  r/(-*7  +  2,n  —  l),n  =  K  —  \,K  —  2, . . .  ,2, 1.  In  the  fifth 
step,  one  can  solve  for  g(J  —  2,  n )  using  the  latent  reference 

point  A  and  the  values  of  r/(J  —  2,  n),n  =  1 . K  -  1.  One 

continues  the  procedure  until  all  the  columns  of g(m.n)  are 
reconstructed,  giving  a  complete  and  unambiguous  recon¬ 
struction  of  g(m,  n ),  and  therefore  of  f(m.  n ). 

If#(l,  K)  ^  0,  then  one  can  alternatively  use  that  point  as 
B  and  perform  the  reconstruction  as  described  above,  but 
reversing  the  roles  of  the  rows  and  columns. 

It  was  recently  noted  that  Eisenstein's  theorem  allows  for 
the  rectangular  region  of  support  (see  Fig.  1 )  to  be  extended 
over  (in  the  same  column  as)  point  .4.  However,  in  that  case, 
there  is  no  simple  recursive  algorithm  for  reconstructing  the 
object. 

B.  Support  Uniqueness  for  Fiddy-Brames-Dainty 
Objects 

In  the  reconstruction  method  described  above,  it  was  im¬ 
plicitly  assumed  that  the  support  of  the  object  function  was 
known.  However,  as  will  be  shown  by  what  follows,  such  an 
assumption  is  not  necessary,  since  an  FBD  object  can  be 
shown  to  be  an  FBD  object  from  its  autocorrelation.  In  order 
to  use  theorems111  relating  to  reconstructing  the  support  of  an 
object  from  the  support  of  its  autocorrelation  function,  during 
the  discussion  of  the  object  and  autocorrelation  supports  we 
assume  that  the  object  function  is  real  and  nonnegative.  (It 
might  happen  that  what  follows  may.  with  appropriate  mod¬ 
ifications,  also  be  true  for  complex-valued  objects:  but  t hi> 
would  require  further  development.) 

Given  only  the  support  of  the  autocorrelation,  one  can 
usually  only  put  an  upper  bound  on  the  support  of  the 
object."’  Such  upper  bounds,  sets  that  can  contain  translates 
of  the  supports  of  all  possible  solutions,  we  refer  to  as  locator 
sets.  One  such  locator  set  is  the  intersection  of  the  autocor¬ 
relation  support  with  a  translate  of  itself,  where  the  translate 
is  such  that  the  center  of  the  second  autocorrelation  support 
is  within  the  first  autocorrelation  support.1"  Assuming  that 
ry(— </  +  1,K  —  1)  *  0.  and  translating  the  one  autocorrelation 
support  so  that  it  is  centered  at  t—J  +  l.K  —  1 ),  one  arrives 
at  the  locator  set  shown  in  Fig.  2  lor  the  case  of  the  FBD  object 
support  shown  in  Fig.  1(b),  In  addition,  since  the  autocor¬ 
relation  is  2-7  +  1  pixels  wide  and  'IK  +  1  pixels  high,  the  object 
must  be  -7  +  1  pixels  wide  and  K  +  1  pixels  high.  Since  the 
object  support  must  be  contained  within  the  locator  set  shown 
in  Fig.  2,  which  is  J  +  2  pixels  wide  and  K  +  2  pixels  high,  the 
object  support  must  include  either  the  lower  left  point  or  t he 
upper  right  point  but  not  both  Keeping  either  one  of  these 
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Fig.  2.  Locator  set  containing  all  possible  solutions,  used  to  show 
that  the  support  solution  is  unique. 


<b) 


Fig  2.  Alternative  case,  (a)  Object  support;  (b)  autocorrelation 
support:  let  locator  set. 

two  points  and  discarding  the  other,  one  is  left  with  the  sup¬ 
port  of  the  object  (or  the  180°  rotated  version— the  twofold 
ambiguity).  Suppose,  on  the  other  hand,  that  r/(-J  +  1,  K 
-  1 )  =  0.  For  example,  suppose  that  the  object  support  is  that 
shown  in  Fig.  2(a).  Then  the  autocorrelation  support  is  that 
shown  in  Fig.  2(b).  A  locator  set,  formed  by  taking  the  in¬ 
tersection  of  this  autocorrelation  support  with  one  translated 
to  be  centered  at  the  first  nonzero  point  in  row  (K  —  1),  is 
shown  in  Fig.  2(c).  As  in  the  case  of  Figs.  1  and  2,  since  the 
autocorrelation  is  IK  +  1  pixels  high,  the  object  must  be  K  + 
1  pixels  high,  and  therefore  either  the  lower  right  or  the  upper 
left  point  (but  not  both)  in  Fig.  3(c)  must  be  within  the  object 
sup|Mirt.  Suppose  we  take  the  lower  left  point  as  being  within 
the  object  (choosing  the  upper  right  point  will  result  in  the 
180°  rotated  solution).  Then,  since  the  autocorrelation  is  2J 
+  1  pixels  wide  and  therefore  the  object  must  be  J  +  1  pixels 
wide,  the  object  must  be  contained  within  the  first  J  +  1  col¬ 
umns  on  the  left  of  Fig.  2(c),  which  is  just  the  support  of  the 
object  as  shown  in  Fig.  2ia).  From  these  examples  it  can  be 
seen  that,  in  general,  if  the  object  is  an  FBI)  object,  then  its 
support  can  be  reconstructed  from  the  autocorrelation  func¬ 
tion.  from  which  it  is  also  evident  that  the  object  has  an  FBI) 
support. 

C.  Triangular  Objects 

Other  tvpes  of  objects,  in  addition  to  FBI)  objects,  can  be  re¬ 
constructed  by  the  recursive  method.  In  this  and  the  next 
section  the  reconstruction  of  two  other  classes  of  objects  are 
shown.  Consider,  tor  example,  objects  whose  support  is 
contained  in  the  triangular  shape  shown  in  Fig.  4(a).  As¬ 
suming  that  the  object's  support  i-  known  a  prmn.  it  has  been 


shown  that  for  this  particular  object  shape  the  boundaries  can 
be  reconstructed  in  a  simple  way,14  assuming  A,  B,  ('  *  0. 
Since  the  vector  spacings  between  points  A  and  B,  B  and  C, 
and  C  and  A  are  each  unique,  from  the  corner  points  in  the 
autocorrelation,  as  shown  in  Fig.  4(b),  we  have 

r(0,  K)  =  /(0,  fO/*(0, 0)  =  C'A*,  (16a) 

rU,  - K )  =  f(J,  0)/*(0,  K)  =  BC *,  ( 16b) 

r(J,  0)  =  f(J ,  0)/*(0,  0)  =  BA*.  (16c) 

Combining  these  gives 

,  r*(0,  K)r{J,  0) 

HJ.-K,  '  ',7' 

Without  loss  of  generality  the  phase  of  A  can  be  chosen  to  be 
zero,  and  then  A  is  given  by  the  positive  square  root  of  Eq. 
(17).  Then  we  can  also  compute 

B  =  r(J,0)/A *,  (18a) 

C  =  r(0,  K)!A*.  (18b) 

Then  the  values  of  the  leftmost  column  of  the  object  are  given 
by 

/(0,n)  =  r(-J,n)/R*.  (19) 

the  values  of  the  bottom  row  are  given  by 

f(m,  0)  =  rim,  -K)/C*,  (20) 

and  the  values  of  the  diagonal  are  given  by 

f(m.  K  -  m)  -  rim,  K  ~  m)/A*.  (21) 

From  this  point  one  could  determine  the  remainder  of  the 
object  by  solving  systems  of  equations, S14  but  an  easier  way 

Kf"*C=l(0,K) 


M  M],-K)  =  BC‘ 


F14:.  I  Triangular  shaped  obted  <.»)  Object  support;  ibi  antm  <>r 
relation  support  The  object  i"  unnjueK  i.imnnj;  1  rintmuiar-'li.iped 
"Mbit ion> i  reconstructed  Iront  its  autororrelat ion  hind  ion. 
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Fig.  S.  Specific  triangular-shaped  object,  (a)  The  object;  (b)  a 
second  nontriangular-shaped  solution;  (c)  the  common  autocorrela¬ 
tion  function;  (d)  the  function  used  to  synthesize  objects  shown  in  (a) 
and  (b). 

is  possible  if  one  cleverly  chooses  the  order  in  which  the 
equations  are  solved.  In  particular,  only  one  linear  equation 
with  one  unknown  at  a  time  need  be  solved,  and  the  solution 
at  each  step  is  unique,  if  one  solves  in  the  following  order.  In 
a  similar  manner  as  was  done  for  the  FBD  objects,  solve  for 
the  points  in  column  m  =  1  using  R  as  a  latent  reference  point, 
and  solve  for  the  points  in  row  n  =  1  using  C  as  a  latent  ref¬ 
erence  point.  Next  solve  for  the  points  in  column  m  =  2  using 
flasa  latent  reference  point,  and  solve  for  the  points  in  row 
n  =  2  using  C  as  a  latent  reference  point.  This  procedure  is 
continued  until  all  of  f(m,  nl  is  reconstructed.  Other  or¬ 
derings  for  the  recursive  solution  of  the  equations  are  also 
possible. 

The  solution  given  above  for  the  triangular-shaped  object 
is  unique  among  objects  having  that  support  but  may  not  be 
unique  among  all  objects.  Momentarily  restricting /(m,  n ) 
to  the  case  of  nonnegative  objects,  one  can  use  the  autocor¬ 
relation  support  tri-intersection  reconstruction  for  convex 
sets1"  to  show  that  there  exists  a  family  of  object  supports  that 
have  autocorrelation  supports  equal  to  the  one  shown  in  Fig. 
4(b).  One  member  of  that  family  is  the  original  object  support 
shown  in  Fig.  4(a).  Another  member  is  an  object  support 
resembling  the  autocorrelation  support  shown  in  Fig.  4(b)  but 
only  half  its  size.  For  these  latter  members  there  is  no  simple 
recursive  reconstruction  algorithm  as  there  is  for  the  trian¬ 
gular-shaped  object. 

Further  insights  can  be  ohtained  by  analyzing  a  simple  case. 
A  case  for  which  there  are  exactly  two  different  solutions  (not 
counting  18()°-rotated  versions)  can  be  obtained  by  starting 
with  nonsymmetric  functions  h\(x,  v)  and  h>(x,  v)  whose 
Fourier  transforms  are  nonfactorable  and  generating  a  first 
object,  which  is  h  t(x.  y )  convolved  with  h  >(x,y),  and  a  second 
object,  which  is  h|(x.y)  convolved  with  h  g(  —  x,  — y)  (i.e.,  the 
cross  correlation). 1  ‘  Two  such  objects,  their  common  auto¬ 
correlation  function,  and  the  h  \(x.y )  =  h  > (x,  v)  used  to  gen¬ 
erate  them  are  shown  in  Figs.  5(a)  through  5(d),  respectively. 
In  this  case  one  obtains  the  "unique”  solution  shown  in  Fig, 
5ta)  it  triangular  support  is  assumed,  and  the  "unique"  solu¬ 
tion  shown  in  Fig  5(b)  if  the  only  other  possible  support  is 
assumed. 

Since  relatively  lew  2  1)  objects  have  factorable  Fourier 
transforms,  the  ambiguous  example  shown  in  Fig.  5  is  unusual. 
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If  one  started  with  a  random  object  having  the  same  support 
as  the  object  in  Fig.  5(b),  and  if  one  incorrectly  assumed  that 
the  object  had  the  same  triangular  support  as  the  object  in  Fig. 
5(a),  then  one  would  obtain  what  at  first  glance  would  appear 
to  be  a  triangular-shaped  solution.  In  the  process  of  calcu¬ 
lating  the  solution  one  would  use  only  the  points  on  the  pe¬ 
rimeter  of  the  autocorrelation  function,  with  which  the  “so¬ 
lution"  would  be  consistent.  However,  on  further  inspection 
one  would  usually  find  that  the  triangular-shaped  solution  is 
inconsistent  with  the  interior  points  of  the  autocorrelation 
function.  Only  in  the  unlikely  event  that  the  original  object’s 
Fourier  transform  is  factorable  would  the  triangular-shaped 
solution  be  completely  consistent  with  the  autocorrelation 
function.  Therefore  if  the  given  autocorrelation  function 
admits  to  a  possible  solution  by  the  recursive  method,  then 
one  should  reconstruct  the  solution  with  the  assumed  support, 
then  compute  its  autocorrelation  function  and  compare  it  with 
the  given  autocorrelation  function  to  determine  whether  the 
assumed  support  is  valid. 

D.  Another  Case 

For  a  final  example,  consider  objects  contained  within  the 
support  shown  in  Fig.  6(a).  Comparing  it  with  Fig.  1(b),  it 
would  be  a  FBD  object  if  it  were  not  for  the  fact  that  R  =  0. 
Assuming  that  the  support  of  the  object  is  known,  it  can  be 
reconstructed  by  the  following  recursive  steps  if  points  A  and 

R'  ^  0  and  if  either  point  C'  or  C"  9*  0.  First  f(J,  2) . f(J , 

K)andf(2,  K), . . .  ,f(J  -  1,  A)  are  solved  using  A  as  the  ref¬ 
erence  point.  A  can  be  determined  from  an  equation  similar 
toEqs.  (11)  and  (17).  Next  C  =  /( 1,  K  -  1),  then  /(l,  K  -  2), 
. . . ,  then  f (1.  2)  are  solved  using  R'  as  the  latent  reference 
point.  Next /(J  -1,1)  is  solved  using  C  or  C"  as  the  latent 
reference  point.  Next  /( 1 ,  1)  is  solved  using  R'  as  the  latent 
reference  point.  Next  f{J  -  1.  2),  ....  f(J  -  1,  K  -  1)  are 
solved  using  A  as  the  latent  reference  point.  Then  the  pattern 

repeats:  solve  for  ^(2  ,K  -  1) . /( 2,  2)  recursively  using  R'. 

then  solve  for  f(J  —  2, 1)  using  C  or  C",  then  solve  for  /( 2,  1 ) 


(a)  (b) 


(c) 


Fig.  li  Another  c.i-e  related  to  FBI)  objects  tai  Object  support; 
ilii  alternative  Mippnrt  reconstruction;  id  autocorrelation  support 
The  nlueol  is  reconstructed  Irotn  its  autocorrelation  lunction.  with 
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using  B\  then  solve  for  /(J  -  2.  2), . .  .  ,f(J  —  2,  K  —  1 )  using 
A,  etc.,  until  all  the  columns  are  solved. 

The  solution  for  this  object  is  unique  among  objects  having 
support  contained  within  the  support  shown  in  Fig.  6(a). 
Howe  ver,  another  support  may  also  be  possible.  In  a  manner 
similar  to  that  used  in  connection  with  Figs.  1-3,  the  possible 
support  solutions  can  be  narrowed  down  to  those  of  Fig.  6(a) 
and  Fig.  6(b),  given  the  autocorrelation  support  shown  in  Fig. 
6(c).  For  the  support  shown  in  Fig.  6(b)  one  can  reconstruct 
the  object  unambiguously  by  solving  a  proper  sequence  of 
equations  using  latent  reference  points  A,  B,  C,  C\  and  D. 
Therefore,  given  the  autocorrelation  function  whose  support 
is  shown  in  Fig.  6(c),  at  most  two  (and  more  probably  only  one) 
solutions  are  possible,  and  each  can  be  reconstructed  using 
a  simple  recursive  algorithm  depending  on  the  support  shown 
in  either  Fig.  6(a)  or  6(b). 

4.  CONCLUSIONS 

A  simple  recursive  algorithm  has  been  devised  for  recon¬ 
structing  an  object  from  its  autocorrelation  function  (or  its 
Fourier  modulus).  It  works  for  several  types  of  sampled 
objects  having  latent  reference  points,  including  those  satis¬ 
fying  the  conditions  described  by  FBD.  The  manner  in  which 
the  algorithm  results  in  a  unique  solution  constitutes  a  proof 
of  uniqueness  for  FBD  objects  (hut  not  necessarily  for  all 
objects  whose  Fourier  transforms  satisfy  Eisenstein’s  theo¬ 
rem).  One  of  the  principal  lessons  learned  here  is  that  the 
detailed  shape  of  the  boundary  of  an  object  plays  a  crucial  role 
in  determining  the  uniqueness  of  the  solution  to  the  phase- 
retrieval  problem. 

One  might  be  able  to  use  this  method  for  continuous  objects 
(as  opposed  to  inherently  sampled  objects)  if  a  dense  enough 
sampling  of  the  autocorrelation  is  available."’ 

Since  the  algorithm  involves  repeatedly  taking  differences 
and  dividing  by  the  values  of  the  latent  reference  points,  it 
may  be  sensitive  to  noise  and  may  require  latent  reference 
points  having  large  values  for  an  accurate  reconstruction. 
(This  may  be  relate.!  to  the  fact  that  a  large  value  of  .4  was 
required  for  a  successful  reconstruction  using  the  iterative 
Fourier  transform  algorithm.1')  Not  all  the  (nonsymmetric) 
points  in  the  autocorrelation  are  used  by  this  algorithm:  im¬ 
proved  accuracy  should  be  expected  if  t  he  reconstruction  al¬ 
gorithm  were  modified  to  use  also  those  additional  points. 
Those  additional  points  may  also  he  used  to  distinguish 
whether  assumptions  about  the  support  of  the  object  (when 
more  than  one  support  solution  is  possible)  are  valid.  For  the 
best  results  one  should  finish  the  reconstruction  by  using  the 
output  of  this  reconstruction  method  as  the  initial  input  to 
the  iterative  Fourier-transform  algorithm."7  which  finds  a 
solution  that  is  most  consistent  with  both  the  measured  data 
and  the  a  priori  constraints. 

The  reconstruction  algorithm  proposed  here  is  applicable 
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to  only  a  relatively  small  number  of  types  of  objects  How 
ever,  the  approach  of  carefully  selecting  the  order  in  which  the 
equations  are  solved  should  be  helpful  in  the  more  general  list 
of  Dallas's  method  by  limiting  the  growth  ot  the  tree  o!  >olu 
tions.7’ 
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MH9.  Holographic  Reconstruction  with  I.atent  Rrfcrcncc 
Points.*  J.  R  FIENUP,  Environmental  Research  Institute  of 
Michigan,  PO  Bax  8618,  Ann  Arbor,  Michigan  JblOT  -  In  I  he  image 
plane  of  a  Fourier-transform  hologram,  one  finds  the  autocorrelation 
of  the  object -plane  distribution  (the  object  plus  a  reference  point  1. 
which  includes  the  autocorrelation  of  the  object,  the  crosscorrelation 
of  the  object  with  the  reference  point  (i.e.,  the  desired  image),  and  the 
conjugate  image.  When  the  reference  point  is  insufficiently  offset 
from  the  object,  then  a  straightforward  reconstruction  is  frustrated 
by  the  overlap  of  the  desired  image  with  the  autocorrelation  term 
One  then  must  solve  the  Fourier  phase  retrieval  problem  or  equal 
ently  reconstruct  the  object  plane  distribution  from  its  autocorre¬ 
lation  function.  For  certain  cases  there  is  a  unique  relal  ionship  be¬ 
tween  (he  Fourier-plane  intensity  and  the  object  plane  distribution, 
even  when  the  reference  point  1-  close  to  the  object  (hence  it  is  only 
a  latent  reference  point  t.  Examples  of  this  are  object  plane  dotri 
butions  whose  Fourier  transforms  satisfy  Eisenstein's  criterion  1  For 
these  and  certain  other  ty  pes  of  object-plane  distributions,  one  can 
digitally  reconstruct  the  object  plane  distribution  from  its  autocor 
relation  function  using  a  recursive  algorithm  This  also  constitutes 
a  proof  of  the  uniqueness  of  phase  retrieval  for  these  types  of 
object  plane  distributions.  The  algorithm  is  similar  to  Dallas'  al 
gorithm  except  that  it  involves  solving  only  one  linear  equation  at  a 
time.  It  has  applications  in  holography,  astronomy,  and  wave-front 
sensing.  (13  min.) 
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Summary .  An  increasing  body  of  theory  indicates  that  the 
phase  retrieval  problem  usually  has  a  unique  solution  for 
2-D  objects.  In  this  paper  experimental  reconstruction 
results  that  support  the  uniqueness  theory  are  shown. 


1  INTRODUCTION 


In  both  optical  and  radio  astronomy,  sometimes  one  can  accu¬ 
rately  obtain  the  modulus  of  the  Fourier  transform  (i.e.,  the  magnitude 
of  the  complex  visibility  function)  of  an  image,  but  not  the  Fourier 
phase.  In  order  to  obtain  an  image  it  then  becomes  necessary  to  re¬ 
trieve  the  Fourier  phase.  Since  the  autocorrelation  function  can  be 
computed  as  the  inverse  Fourier  transform  of  the  squared  Fourier  modu¬ 
lus,  the  problem  is  equivalent  to  reconstructing  an  image  from  its  auto¬ 
correlation. 

In  this  paper  we  are  concerned  with  the  phase  retrieval  problem  in 
optical  astronomy,  in  which  case  one  cannot  rely  on  such  aids  as  closure 
phase  (Jennison  1958 ).,  However  the  results  shown  here  do  have  relevance 
to  radio  astronomy  as  well. 

Several  methods  have  been  put  forward  for  solving  the  phase  retrieval 
problem  (Liu  &  Lohmann  1973;  Napier  &  Bates  1974;  Frieden  &  Currie  1976; 
Baldwin  &  Warner  1978;  Fienup  1978,  1979,  1982;  Bates  et  al.  1982a).  In 
addition  there  are  a  number  of  reconstruction  techniques  that  depend  on 
the  specific  method  of  data  collection,  for  example,  astronomical 
speckle  interferometry  (Bates  1982b).  Of  the  methods  that  would  work 
for  the  most  general  case,  the  iterative  input-output  Fourier  transform 
algorithm  (Fienup  1978,  1979,  1982)  appears  to  be  the  most  practical. 

when  any  of  the  reconstruction  algorithms  finds  a  solution,  the  question 
remains:  is  it  the  only  (unique)  solution  or  is  it  one  of  many  possible 
(ambiguous)  solutions?  In  Section  2  the  theory  of  the  uniqueness  will 
be  briefly  reviewed.  Then  in  Section  3  experimental  reconstruction 
results  will  be  shown  that  are  consistent  with  the  theory  that  the  2-D 
case  is  usually  unique.  In  addition,  experimental  reconstruction 
results  are  shown  that  indicate  that  noise  in  the  Fourier  modulus  data 
does  not  radically  change  the  uniqueness  of  the  solution. 
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2  UNIQUENESS  THEORY 

When  we  speak  of  the  reconstruction  being  unique  or  ambig¬ 
uous,  we  ignore  translations  and  180°  rotations  since  neither  of  these 
operations  affects  the  Fourier  modulus.  Here  we  are  also  assuming  that 
the  object  has  a  finite  spatial  (or  angular)  extent. 

The  one-dimensional  ( 1 — D )  phase  retrieval  problem  has  long  been  known  to 
be  highly  ambiguous  (Walther  1963).  Only  for  the  special  cases  of 
objects  known  to  consist  of  sufficiently  separated  parts  or  nonnegative 
objects  having  sufficiently  separated  parts  is  the  1-D  phase  retrieval 
problem  usually  unique  (Greenaway  1977;  Crimmins  &  Fienup  1983). 

The  2-D  case  is  quite  different.  This  can  best  be  understood  from  the 
theory  developed  by  Bruck  and  Sodin  (1979).  They  considered  the  special 
case  of  an  object  sampled  on  a  rectangular  lattice.  For  the  1-D  case 
the  Fourier  transform  can  then  be  expressed  as  a  polynomial  of  order  M 
of  a  single  complex  variable,  and  such  a  polynomial  can  always  be  fact¬ 
ored  into  M  irreducible  factors  (by  the  fundamental  theorem  of^a^gebra). 
They  showed  that  this  implies  that  in  the  1-D  case  there  are  2  possi¬ 
ble  solutions,  although  not  all  of  those  solutions  would  satisfy  a  non¬ 
negativity  constraint  (Bates  1969).  On  the  other  hand,  polynomials  of 
two  complex  variables  having  arbitrary  coefficients  are  only  rarely 
factorable.  Consequently  the  2-D  case  is  usually  unique.  Although  the 
2-D  theory  for  continuous  functions  has  not  yet  been  fully  developed,  it 
is  likely  that  a  similar  result  will  hold. 

Of  course  one  can  always  fabricate  2-D  examples  that  are  not  unique.  An 
example  is  an  object  formed  by  convolving  two  nonnegative  functions.  A 
second  object,  formed  by  convolving  the  first  nonnegative  function  with 
an  inverted  (i.e.,  rotated  by  180°)  version  of  the  second  nonnegative 
function,  has  the  same  Fourier  modulus  as  the  first  object.  Another 
method  of  synthesizing  ambiguous  cases  was  given  by  Huiser  and  van  Toorn 
(1980).  However,  these  fabricated  ambiguous  objects  are  very  special 
cases--most  2-D  objects  do  not  fit  into  these  categories. 

There  are  also  a  number  of  classes  of  objects  for  which  the  phase  re¬ 
trieval  problem  is  known  to  be  unique  (as  opposed  to  just  being  usually 
unique).  For  example,  if  the  object  includes  an  unresolved  (delta- 
function-like)  point  far  enough  away  from  the  rest  of  the  object,  then 
the  autocorrelation  includes  the  rest  of  the  object  as  one  of  its  terms 
(Liu  &  Lohmann  1973),  analogous  to  holography.  It  has  also  been 
recently  discovered  that  for  objects  having  a  special  support  there  is  a 
unique  reconstruction  even  if  the  reference  points  are  very  close  to  the 
rest  of  the  object  (Fiddy  et  al.  1983).  The  support  of  an  object  is  the 
set  of  points  over  which  it  is  nonzero,  i.e.,  its  shape.  Also  using 
latent  reference  points  it  can  be  shown  that  these  and  other  objects 
having  certain  supports  can  be  uniquely  reconstructed  from  their  Fourier 
modulus  (Fienup  1983a).  These  recent  results  point  to  the  importance  of 
the  support  of  an  object  in  determining  whether  the  object  can  be 
uniquely  reconstructed  from  its  Fourier  modulus.  Methods  for  recon- 
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structing  support  information  without  resorting  to  a  complete 
reconstruction  are  also  being  investigated  (Fienup  et  al.  1982). 

3  UNIQUENESS  EXPERIMENTS 

3 . 1  Iterative  reconstruction  algorithm 


One  approach  to  determining  whether  most  objects  of  interest 
are  uniquely  reconstructable  from  their  Fourier  modulus  is  to  perform  a 
number  of  reconstruction  experiments.  This  is  now  possible  due  to  the 
existence  of  a  practical  reconstruction  algorithm,  the  iterative  Fourier 
transform  algorithm  (Fienup  1978,  1979,  1982). 

The  iterative  Fourier  transform  algorithm  uses  all  the  available  mea¬ 
surements  and  a  priori  information  to  arrive  at  a  solution.  In  the 
Fourier  domain  one  has  the  measured  Fourier  modulus  data,  which  is  an 
estimate  of  the  true  modulus  of  the  Fourier  transform  of  the  object.  In 
the  object  domain  one  has  the  a  priori  constraint  that  the  object's 
spatial  (or  angular)  brightness  distribution  is  a  nonnegative  function. 
From  the  Fourier  modulus  data  one  can  compute  an  estimate  of  the  ob¬ 
ject's  autocorrelation  function.  From  the  autocorrelation  one  can  place 
upper  bounds  on  the  diameter  of  the  object  (only  in  special  cases  can 
the  support  of  the  object  be  readily  determined  from  the  support  of  its 
autocorrelation)  (Fienup  et  al.  1982). 

The  iterative  Fourier  transform  algorithm  is  a  modification  of  the 
Gerchberg-Saxton  (1972)  algorithm  that  has  been  used  in  electron  micros¬ 
copy  and  for  other  applications  (Fienup  1983b).  The  simplest  version  of 
the  iterative  algorithm  consists  of  the  four  following  steps.  (1)  An 
estimate  of  the  object  (an  input  image)  is  Fourier  transformed.  (2)  The 
resulting  Fourier-domain  function  is  forced  to  conform  to  the  measure¬ 
ments  by  replacing  the  computed  Fourier  modulus  with  the  measured  Four¬ 
ier  modulus.  (3)  The  result  is  inverse  Fourier  transformed,  yielding  an 
output  image.  (4)  A  new  input  image  is  formed  by  forcing  the  output 
image  to  conform  to  the  object-domain  constraints,  i.e.,  it  is  set  equal 
to  zero  where  it  is  negative  or  where  is  exceeds  the  known  diameter 
(i.e.,  the  support  constraint).  This  algorithm,  which  we  call  the 
error-reduction  algorithm,  can  be  proven  to  converge  in  the  sense  that 
the  error  at  t|ig  k  iteration  is  always  less  than  or  equal  to  the  error 
at  the  (k  -  1)  iteration.  Here  the  error  is  defined  as  the  amount  by 
which  the  computed  Fourier  modulus  differs  from  the  measured  Fourier 
modulus  or  as  the  amount  by  which  the  output  image  violates  the 
object-domain  constraints.  However,  in  practice  the  error-reduction 
algorithm  usually  converges  so  slowly  <_■  at  it  is  impractical  for  this 
application  (Fienup  1982). 

Fortunately  there  exist  a  number  of  accelerated  versions  of  the  algo¬ 
rithm  which  converge  in  a  reasonable  number  of  iterations.  To  date  the 
fastest  version  of  the  algorithm  is  the  hybrid  input-output  algorithm. 
Its  first  three  steps  are  identical  to  those  of  the  error-reduction 
algorithm  described  above.  The  fourth  step  of  the  hybrid  input-output 
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algorithm  consists  of  forming  a  new  input  image  that  is  equal  to  the 
output  image  wherever  the  output  image  satisfies  the  constraints,  and  is 
equal  to  the  previous  input  image  minus  a  constant  factor  times  the 
output  image  wherever  the  output  image  violates  the  constraints.  Any 
value  between  0.5  and  1.0  works  well  for  constant  factor,  which  is 
similar  to  a  negative  feedback  parameter. 

In  one  series  of  trials,  the  algorithm  was  run  on  a  fabricated  Fourier 
modulus  which  was  known  to  have  two  solutions.  One  of  the  two  solutions 
was  reconstructed  in  about  half  of  the  trials  and  the  other  solution  was 
reconstructed  in  the  other  half  of  the  trials.  Which  of  the  two  solu¬ 
tions  was  obtained  depended  on  the  array  of  random  numbers  used  as  the 
initial  input  to  the  algorithm.  Therefore  we  believe  that  if  there  are 
multiple  solutions,  then  the  algorithm  is  equally  likely  to  find  any  one 
of  them  (if  the  initial  input  is  sufficiently  random  and  unbiased),  and 
if  run  enough  times  with  different  initial  inputs,  it  will  probably  find 
all  of  them.  In  a  practical  reconstruction  situation  in  which  the  solu¬ 
tion  is  not  known  beforehand,  if  one  were  to  run  the  algorithm  two  or 
three  times,  each  time  using  a  different  array  of  random  numbers  for  the 
initial  input,  and  if  the  reconstructed  images  were  the  same  each  time, 
then  one  would  be  highly  confident  that  one  had  found  the  solution  and 
that  it  is  unique  (Fienup  1979). 

A  problem  with  experimental  reconstruction  experiments  is  that  there  is 
no  guarantee  that  the  iterative  algorithm  will  converge  to  any  solution, 
even  when  an  accelerated  version  of  the  algorithm  is  used.  One  can 
ttjink  of  the  reconstruction  algorithm  as  an  iterative  search  through  an 
N  -dimensional  parameter  space  (each  dimension  or  parameter  correspond¬ 
ing  to  the  value  of  one  of  the  pixels  of  the  image),  seeking  to  minimize 
the  error  of  the  estimate.  While  searching  for  the  global  minimum  of 
the  error^  the  algorithm  could  stagnate  at  a  local  minimum  of  the  error 
in  that  N  -dimensional  space.  The  likelihood  of  stagnation  and  the 
success  of  the  algorithm  depend  on  the  N  -dimensional  topography  of  the 
error  function,  which  varies  from  one  type  of  object  to  another. 
Therefore,  for  particularly  difficult  objects,  i.e.,  ones  for  which  the 
error  has  many  local  minima,  one  may  not  be  able  to  test  for  uniqueness 
since  the  reconstruction  algorithm  fails.  Fortunately  such  a  problem 
has  occurred  only  occasionally  for  the  types  of  objects  examined. 

One  particular  convergence  problem  has  occurred  on  several  occasions. 
Sometimes  the  algorithm  stagnates  at  a  deep  local  minimum  at  which  the 
output  image  resembles  the  original  object  but  with  a  pattern  of  stripes 
superimposed.  A  similar  phenomenon  has  occurred  in  other  reconstruction 
situations  (Cornwell  1983).  In  most  cases  the  stripes  are  of  low  con¬ 
trast,  superimposed  on  an  otherwise  excellent  reconstructed  image,  and 
are  of  little  concern.  In  other  cases  the  stripes  are  of  high  enough 
contrast  to  be  objectionable.  When  the  Fourier  modulus  data  is  suffi¬ 
ciently  noisy,  then  the  stripes  do  not  appear  (Feldkamp  &  Fienup  1980). 
The  nature  of  the  stripes  is  as  yet  not  fully  understood  and  methods  of 
avoiding  them  remain  to  be  developed. 
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An  example  of  the  stripes  phenomenon  is  shown  in  Figure  1.  Figure  1(a) 
shows  the  original  object  and  Figure  Kb)  shows  a  reconstructed  image, 
which  appears  to  be  quite  faithful.  Figure  1(c)  shows  the  same  recon¬ 
structed  image,  but  heavily  overexposed,  in  order  to  emphasize  the 
low-contrast  vertical  stripes  that  are  present,  although  difficult  to 
discern,  in  the  image.  Figures  l(d-f)  show  the  overexposed  reconstruct¬ 
ed  images  resulting  from  three  other  trials  of  the  algorithm,  each  of 
which  was  initialized  with  a  different  array  of  random  numbers.  In 
each  of  these  three  cases  the  reconstructed  image  contains  a  more  easily 
discernable  pattern  of  stripes,  but  the  spatial  frequencies  and  orienta¬ 
tions  of  the  stripes  are  different  in  each  case.  The  stripes  extend 
throughout  image  space  (although  they  are  weaker  away  from  the  support 
of  the  object),  and  therefore  by  inspection  of  the  reconstructed  images 
it  is  possible  to  determine  that  the  stripes  are  an  artifact  rather  than 
a  true  feature  of  the  object.  Furthermore,  it  is  possible  to  discern 
the  true  image  from  the  stripes  since  the  stripes  change  from  one  recon¬ 
struction  to  the  next,  but  the  true  features  of  the  object  are  present 
in  all  the  reconstructed  images. 

3.2  Experimental  uniqueness  results  for  various  objects 

The  iterative  reconstruction  algorithm  was  used  to  recon¬ 
struct  a  number  of  different  objects  from  their  Fourier  modulus.  The 
objects  examined  are  of  a  very  practical  and  interesting  class: 
digitized  photographs  of  satellites.  They  also  share  a  feature  that  we 
suspect  makes  them  "good"  objects  to  reconstruct:  they  have  interesting 
(i.e.,  complicated)  shapes. 

A  typical  result  is  shown  in  Figure  2,  in  which  (a)  is  the  original 
object  and  (b)  is  the  reconstructed  image  (Fienup  1981).  For  this  and 
almost  all  of  the  cases  examined,  the  reconstructed  image  looks  much 
like  the  original  object  except  for  differences  that  could  be  attributed 
to  stripes.  For  example,  horizontal  stripes  are  evident  over  portions 
of  the  reconstructed  image  shown  in  Figure  2(b).  Therefore,  except  for 
the  presence  of  the  stripes  artifact  which  we  believe  is  a  character¬ 
istic  of  a  local  minimum  rather  than  an  inherent  ambiguity,  most  objects 
of  this  type  are  uniquely  related  to  their  Fourier  modulus. 

There  are  exceptions,  however.  Figure  3  shows  one  case  that  worked 
particularly  poorly.  The  object  shown  in  Figure  3(a)  is  nearly  centro- 
symmetric.  Figure  3(b)  shows  the  reconstructed  image,  which  is  not  very 
faithful.  This  particular  case  has  similarities  with  the  ambiguous  case 
fabricated  by  Huiser  and  van  Toorn  (1980).  From  this  we  see  that,  al¬ 
though  ambiguous  cases  may  be  unusual,  they  are  by  no  means  nonexistent 
in  the  real  world. 

3 . 3  Experimental  uniqueness  in  the  presence  of  noise 

As  with  any  reconstruction  method  the  sensitivity  of  the 
algorithm  to  noise  is  a  major  point  of  concern.  Reconstruction  results 
using  noisy  Fourier  modulus  data  have  shown  that  the  iterative  Fourier 
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Figure  1.  (a)  Original  object;  (b)  image  reconstructed 

from  Fourier  modulus  using  iterative  algorithm;  (c)-(f) 
four  images  reconstructed  using  different  starting 
inputs — these  pictures  were  intentionally  overexposed  in 
order  to  emphasize  the  stripes. 
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Figure  2.  la)  A  typical 
object;  (b)  image  recon¬ 
structed  from  Fourier 
modulus . 


Figure  3.  ( a )  An  atypical 

object,  for  which  the  re¬ 
constructed  image  (b)  does 
not  resemble  the  object. 
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transform  algorithm  is  not  highly  sensitive  to  noise  (Fienup  1978, 

1979).  In  this  section  the  results  of  a  systematic  study  of  the  noise 
sensitivity  of  the  reconstruction  (Feldkamp  &  Fienup  1980)  are 
summarized. 

When  noise  is  present  in  the  Fourier  modulus  data,  then  there  is 
generally  no  solution  that  is  completely  consistent  with  both  the 
measured  data  and  the  constraints.  For  example,  an  autocorrelation 
function  computed  from  a  noisy  Fourier  modulus  would  be  very  likely  to 
have  some  negative  values  for  the  largest  separations.  Obviously  no 
nonnegative  object  can  have  an  autocorrelation  having  negatives;  there¬ 
fore  there  could  be  no  nonnegative  object  consistent  with  the  noisy 
Fourier  modulus.  Nevertheless  the  algorithm  searches  for  a  solution 
that  is  most  consistent  with  the  measured  data  and  constraints,  and  in 
doing  so  it  can  arrive  at  a  useful  image. 

Fourier  modulus  data  was  simulated  to  have  the  type  of  noise  that  would 
be  present  in  astronomical  speckle  interferometry.  The  object  shown  in 
Figure  1(a)  was  convolved  with  156  different  point-spread  functions  to 
produce  156  different  blurred  images.  Each  of  the  point-spread  func¬ 
tions  represents  a  different  realization  of  the  blurring  due  to  the 
turbulent  atmosphere.  The  widths  of  the  point-spread  functions  were 
comparable  to  the  width  of  the  object.  The  blurred  images  were  then 
subjected  to  a  Poisson  noise  process  to  simulate  the  effects  of  photon 
noise.  The  degraded  images  were  then  processed  to  produce  a  noisy 
Fourier  modulus  estimate  by  Labeyrie's  (1970)  method,  as  modified  by 
Goodman  and  Belsher  (1976)  to  eliminate  the  bias  noise  term  from  the 
squared  Fourier  modulus. 

Figure  4  shows  a  noise-free  Fourier  modulus  (a)  and  three  examples  of 
the  simulated  noisy  Fourier  modulus  estimates  (b)-(d)  with  increasing 
noise.  Figure  5  shows  the  original  undegraded  object  (a)  and  three 
images  (b)-(d)  reconstructed  from  the  respective  noisy  Fourier  modulus 
estimates  of  Figure  4.  For  the  case  shown  in  Figures  4(b)  and  5(b), 
which  represent  a  realistic  amount  of  noise  for  this  situation,  the 
normalized  rms  error  of  the  Fourier  modulus  estimate  was  2.9%  and  the 
reconstructed  image  is  very  good.  For  the  case  shown  in  Figures  4(c) 
and  5(c),  only  1/50  as  many  photons  were  assumed  to  be  available,  and 
the  rms  error  of  the  Fourier  modulus  estimate  is  a  very  poor  32%; 
nevertheless  the  reconstructed  image  still  retains  some  recognizable 
features.  In  the  case  shown  in  Figures  4(d)  and  5(d),  an  extreme  amount 
of  noise  was  present,  and  the  rms  error  of  the  Fourier  modulus  estimate 
is  near  100%;  since  this  Fourier  modulus  estimate  does  not  resemble  the 
true  Fourier  modulus,  then,  as  one  would  expect,  the  reconstructed  image 
does  not  resemble  the  original  object. 

4.  CONCLUSIONS 


Theory,  which  points  toward  the  conclusion  that  a  2-D  object 
of  finite  extent  is  ordinarily  uniquely  related  to  the  modulus  of  its 
Fourier  transform,  has  been  supported  by  experimental  reconstruction  re- 
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Figure  4.  Fourier  modulus  estimates  with  noise,  having  rms 
error  (a)  0%,  <b)  2.9%,  <c)  32%,  (d)  '100%. 
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suits.  The  vast  majority  of  reconstructed  images  of  satellites  resemble 
the  original  objects  from  which  the  Fourier  modulus  was  computed.  Fur¬ 
thermore,  contrary  to  some  predictions  (Huiser  &  van  Toorn  1980),  the 
uniqueness  properties  do  not  change  radically  when  noise  is  present. 
Rather,  as  more  noise  is  introduced  into  the  Fourier  modulus  estimate, 
the  reconstructed  image  simply  becomes  correspondingly  noisier,  and 
degrades  in  a  gradual  manner. 
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Comments  on  “The  Reconstruction  of  a  Multidimensional 
Sequence  from  the  Phase  or  Magnitude  of  Its 
Fourier  Transform” 
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Abstract-Ylben  one  imposes  a  nonnegativity  constraint,  one  usually 
can  reconstruct  a  two-dimensional  sequence  of  finite  support  from  the 
modulus  of  its  Fourier  transform  using  an  iterative  algorithm,  even 
when  the  initial  estimate  is  an  array  of  random  numbers. 

In  a  recent  paper,1  the  description  of  an  iterative  algorithm 
for  reconstructing  a  sequence  from  the  magnitude  of  its 
Fourier  transform  unintentionally  gives  the  appearance  of  dis¬ 
cussing  an  algorithm  published  earlier  |l|.  In  the  following. 
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the  differences  between  the  algorithm  and  experiments  de¬ 
scribed  by  Hayes1  and  those  published  earlier  [  1  ]  are  clarified. 

Hayes1  reviews  both  the  problem  of  reconstructing  a  se¬ 
quence  from  the  phase  of  its  Fourier  transform  and  the  prob¬ 
lem  of  reconstructing  a  sequence  from  the  magnitude  of  its 
Fourier  transform.  For  the  latter  problem,  he  describes  an 
iterative  algorithm  for  solving  the  problem  as  follows.  “Spe¬ 
cifically,  this  algorithm  involves  the  repeated  Fourier  transfor¬ 
mation  between  the  time  and  frequency  domains  where,  in 
each  domain,  the  known  information  about  the  desired  se¬ 
quence  is  imposed  on  the  current  estimate.  In  the  time  do¬ 
main,. for  example,  a  sequence  is  constrained  to  have  a  given 
region  of  support  whereas  in  the  frequency  domain,  the 
sequence  is  constrained  to  have  a  given  transform  magni¬ 
tude.”1  He  then  shows  examples  where  the  algorithm  de¬ 
scribed  above  fails.  This  failure  should  not  reflect  poorly  on 
the  earlier  work  [  1  ]  since  the  algorithm  described  in  the  quo¬ 
tation  above  and  the  experiments  performed  by  Hayes  differ 
in  important  ways  from  the  earlier  work.  In  Hayes  experi¬ 
ments,  both  the  type  of  information  which  was  assumed  to  be 
known  and  the  reconstruction  algorithm  which  was  used 
differed  from  those  of  the  earlier  work  1 1  ). 

Hayes  is  correct  in  stating1  that  the  magnitude  of  the  Fou¬ 
rier  transform  is  insufficient  to  uniquely  specify  a  sequence; 
additional  information  or  constraints  are  required.  Depending 
on  the  application,  one  often  has  available  additional  informa¬ 
tion  or  constraints,  and  a  reconstruction  may  then  be  possible 
[  2),  ( 3],  Two  important  constraints  which  often  occur  (as  in 
astronomy)  are  a  known  support  (or  bounds  on  the  support) 
of  a  sequence,  and  the  constraint  that  the  sequence  be  nonneg¬ 
ative  [4).  Unlike  the  algorithm  used  by  Hayes,1  the  iterative 
method  described  earlier  [  1  ]  primarily  uses  the  nonnegativity 
constraint.  Using  the  iterative  algorithm,  we  have  been  very 
successful  in  reconstructing  two-dimensional  nonnegative  se¬ 
quences  from  their  Fourier  magnitude  [l]-[6]  In  this  case, 
the  sequences  must  have  finite  support,  but  it  is  possible  to 
reconstruct  them  even  when  the  support  is  not  known.  Ex¬ 
cept  for  special  cases,  it  is  not  possible  to  determine  the  sup¬ 
port  of  a  sequence  from  the  support  of  its  autocorrelation 
(which  is  the  inverse  Fourier  transform  of  the  squared  Fourier 
magnitude)  [7],  so  the  support  information  is  usually  not 
available  anyway  One  can  only  place  upper  bounds  on  the 
support  (7).  If  an  upper  bound  on  the  support  is  utilized  dur¬ 
ing  the  iterations,  then  the  algorithm  converges  faster  (in 
about  100  or  200  iterations  for  our  work)  than  when  using 
only  the  nonnegativity  constraint  (in  which  case  we  found  that 
several  hundred  iterations  are  required) 

Unlike  the  algorithm  used  by  Hayes,  the  iterative  algorithm 
described  earlier  [  1  |  does  not  simply  satisfy  the  constraints 
(nonnegativity  and  bounds  on  the  support)  in  the  time-domain 
step  of  the  iteration  Such  an  algorithm,  which  we  refer  to  as 
the  error-reduction  algorithm,  was  discussed  earlier  [  1  |  where 
it  is  noted  that.  "For  the  present  application,  the  error-reduc¬ 
tion  approach  requires  an  impractically  large  number  of  itera¬ 
tions  for  convergence."  It  is  only  a  version  of  the  input-output 
algorithm  (  1  |-[6|  which  is  capable  of  converging  in  100  or  so 
iterations. 

Hayes  found  that  if  the  initial  estimate  used  in  the 
iteration  has  a  Fourier  transform  with  the  correct  magnitude 
and  either  zero  phase  or  random  phase,  then  the  iteration  will 
not  generally  converge  to  the  correct  sequence.”1  However, 
using  the  input-output  algorithm  with  a  nonnegativity  con¬ 
straint.  we  obtained  good  reconstruction  results  when  the  al¬ 
gorithm  was  initialized  with  arrays  of  random  numbers  (  1  |- 
[6 1  The  algorithm  has  also  been  shown  to  be  surprisingly 
insensitive  to  noise  [5], 

When  the  error-reduction  algorithm  was  used  with  a  non¬ 


negativity  constraint  (as  well  as  a  support  constraint),  it  took 
many  thousands  of  iterations  for  convergence  |3|.  (6|.  There 
fore,  if  one  were  to  employ  the  error-reduction  algorithm 
without  a  nonnegativity  constraint,  then  one  would  expect 
convergence  to  take  much  longer,  if  it  ever  converges.  Conse¬ 
quently,  it  is  consistent  with  our  experience  that  the  type  of 
reconstruction  experiments  performed  by  Hayes  would  be 
unsuccessful. 

Of  course,  there  are  situations  for  which  the  nonnegativity 
constraint  does  not  apply.  Then  one  might  wonder  whether 
it  is  possible  to  reconstruct  a  sequence  of  finite  support  from 
its  Fourier  magnitude.  Theory  ([81,  Hayes1  )  seems  to  indicate 
that  the  solution  will  usually  be  unique.  However,  as  shown 
by  Hayes,  the  error-reduction  algorithm  is  not  a  practical  ap¬ 
proach  to  finding  the  solution.  One  might  possibly  succeed 
using  an  accelerated  algorithm,  such  as  the  input-output  al¬ 
gorithm  or  a  gradient  search  method  [61.  but  this  is  an  area 
that  needs  further  work. 

It  should  also  be  noted  that  in  the  phase  retrieval  problem  of 
X-ray  crystallography,  one  reconstructs  the  three-dimensional 
electron  density  function  from  its  Fourier  magnitude  For 
that  problem,  one  has  the  constraints  that  the  electron  density 
is  nonnegative  and  that  it  consists  of  a  discrete  number  of 
atoms.  For  that  problem,  a  number  of  reconstruction  meth¬ 
ods  have  been  developed  [9],  For  the  phase  retrieval  problem 
in  electron  microscopy,  for  which  both  the  wave  function  and 
its  Fourier  transform  are  complex  valued,  one  has  the  addi¬ 
tional  constraint  of  knowing  the  magnitude  of  the  wave  (unc¬ 
tion.  For  that  problem,  the  error-reduction  algorithm  has 
been  shown  to  perform  very  well  [10],  [11] 

In  conclusion,  Hayes’  remark  that  .  .  even  for  those  se¬ 
quences  which  are  uniquely  defined  by  their  magnitude,  it  ap¬ 
pears  that  a  practical  algorithm  is  yet  to  be  developed  for  re¬ 
constructing  a  sequence  from  only  its  magnitude”1  is  stristly 
true  when  no  other  information  is  available;  however,  tor  a 
number  of  important  applications,  there  is  auxiliary  infor¬ 
mation,  such  as  a  nonnegativity  constraint,  and  practical  re¬ 
construction  algorithms  do  exist 
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ABSTRACT 

It  is  shown  that  knowledge  of  the  edges  of  an  object  is  not  always 
sufficient  to  uniquely  reconstruct  an  object  fr^m  the  modulus  of  its 
Fourier  transform  via  the  autocorrelation  function.  On  the  other  hand, 
I  in  some  cases  not  only  can  the  boundary  values  be  determined  from  the 

autocorrelation,  but  also  the  object  can  be  reconstructed  uniquely,  even 
for  complex- valued  objects. 

I 


1.  Introduction 

In  a  number  of  disciplines,  including  astronomy,  x-ray 
crystal  1 ography ,  electron  microscopy,  and  wavefront  sensing,  one 
encounters  the  phase  retrieval  problem.  One  wishes  to  reconstruct 
f(m,  n),  an  object  function,  from  |F(p,  q)|  ,  the  modulus  of  its  Fourier 
transform,  where 

F(p,  q)  =  | F ( p ,  q)  |  exp  [itHp,  q)]  =  n)] 

P-1  Q-l 

■II  f(m,  n)  exp  [-i 2rr(mp/P  +  nq/Q)],  (1) 

m=0  n=0 

where  m,  p  -  0,  1,  . . . ,  P  -  1  and  n,  q  =  0,  1,  ....  Q  -  1.  The  discrete 
transform  is  employed  here  since  in  practice  one  deals  with  sampled  data 
in  a  computer.  The  problem  of  reconstructing  the  object  from  its 
Fourier  modulus  is  equivalent  to  reconstructing  the  Fourier  phase, 

Hp,  q),  from  the  Fourier  modulus,  since  once  one  has  the  phase  as  well 
as  the  modulus,  one  can  easily  compute  f(m,  n)  by  the  inverse  (discrete) 
Fourier  transform.  rf(m,  n),  the  (aperiodic)  autocorrelation  of 
f(m,  n),  is  given  by 

M-l  N-l 

rf(m,  n)  =  ^  £  f(j,  k)f*(j  -  m,  k  -  n)  (2) 

j=0  k=0 

M-l  N-l 

=  X!  X f  (j*  +  m» k  +  n)  (3) 

j=0  k=0 


-  9r_1[iF(p,  q)i2]. 


(4) 


where  the  asterisk  denotes  complex  conjugate  and  where  it  is  assumed 
that  f(j,  k)  =  0  for  m  outside  of  [0,  M  -  1]  and  for  n  outside  of  [0,  N 
-  1].  Note  that  in  order  to  avoid  aliasing  in  the  computation  of 
! F( p »  q )!  it  is  necessary  to  have  M<  P/2  and  N  <  Q/2.  Since  the 
autocorrelation  function  is  easily  computed  from  the  Fourier  modulus  by 
Eq.  (4),  the  phase  retrieval  problem  is  equivalent  to  reconstructing  an 
object  from  its  autocorrelation  function. 

Several  phase  retrieval  algorithms  have  been  proposed,  all  of  them 
requiring  some  additional  measurements  or  constraints  on  the  solution. 

Examples  include  a  reference  point  at  least  one  object-diameter  from  the 

2  3 

object  (giving  rise  to  the  holography  condition  ),  a  second  intensity 

4-5 

measurement  in  another  plane  (in  electron  microscopy  or  wavefront 

6-8 

sensing),  nonnegativity  and  limited  spatial  extent  (in  astronomy), 

Q 

atomic  models  (in  x-ray  crystallography),  objects  consisting  of 
collections  of  points  having  nonredundant  spacings^,  and  objects  having 
latent  reference  points**  (not  satisfying  the  holography  condition). 
For  each  of  these  situations  there  is  a  proof  of  uniqueness  of  the 
solution  that  relies  on  the  types  of  measurements  made,  on  the  a  priori 
information  available,  or  on  the  nature  of  the  reconstruction  algorithm 
i  tsel f . 

Antoher  proposed  phase  retrieval  algorithm  is  a  recursive  one  that 

relies  on  a  priori  knowledge  of  the  boundary  conditions  (i.e.  the  values 

of  the  edges  of  the  object).  The  purpose  of  this  paper  is  to  show 

that  the  general  uniqueness  claims  made  concerning  phase  retrieval  using 

12 

boundary  conditions  are  incorrect;  but  by  the  approach  of  using  latent 
reference  points,**  special  classes  of  objects  can  be  shown  to  be 
unique,  for  complex-valued  objects  as  well  as  for  real-valued  objects. 
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uity  Using  Boundary  Conditions 


2.  Ambiq 

In  Reference  12  a  recursive  algorithm  was  put  forward  for 
reconstructing  an  object  from  the  modulus  of  its  Fourier  transform,  via 
the  autocorrelation  function,  using  boundary  conditions,  i.e.,  assuming 
knowledge  of  the  edges  of  the  object.  A  real-valued  object,  f(m,  n), 
was  assumed  to  be  zero  outside  of  the  rectangular  region  of  support  0< 
ml  M  -  1  and  0  1  n  1  N  -  1.  The  top  and  bottom  nonzero  rows,  3(m)  = 
f(m,  N  -  1)  and  a(m)  =  f(m,  0),  respecti vely ,  and  the  leftmost  and 
rightmost  nonzero  columns,  f(0,  n)  and  f(M  -  1,  n),  respectively,  are 
assumed  to  be  known  a  priori.  Rows  1  and  N  -  2  can  then  be  determined 
by  solving  a  system  of  2M  -  1  linear  equations  in  2M  -  4  unknowns.  For 
example,  from  Eq.  (3)  we  have,  for  n  =  N  -  2,  the  second  from  the  top 
row  of  the  autocorrelation: 


M-l  N-l 

r(m,  N  -  2)  =  II  f  (j>  k)f ( j  +  m,  k  +  N  -  2) 

j=0  k=0 

M-l  M-l 

=  X  f*(j’  0)f(j  +  m,  N  -  2)  +  X  Mj.  +  m,  N  -  1) 

j=0  j=0 

M-l  M-l 

=  X  +  m,  N  -  2)  +  X  Mj,  DB(j  +  m)  (5) 

j=0  j=0 

for  m  =  -M  +  1,  ...,  M-l.  These  are  2M  -  1  equations,  one  for  each 
value  of  m,  in  2M  -  4  unknowns,  f(j,  N  -  2)  and  f(j,  1),  for  j  =  1,  2, 
...,  M  -  2.  Recall  that  a(j),  B(j),  f(0,  N  -  2),  f(M  -  1,  N  -  2), 
f(0,  1)  and  f(M  -  1,  1)  are  assumed  known.  After  f(j,  N  -  2)  and  f ( j , 
1)  are  determined  by  solving  the  system  of  equations  given  in  Eq.  (5) 
above,  then  one  can  solve  for  f(j,  N  -  3)  and  f ( j ,  2)  using  r(m,  N  -  3) 
in  a  similar  manner.  The  remaining  rows  of  the  object  are  solved 


recursively  in  a  similar  manner. 


The  method  described  above  could  work  if  the  systems  of  equations 

have  a  unique  solution  for  the  unknowns.  Restricting  the  solution  to 

12 

real-valued  f  s,  a  claim  has  been  made  that,  "It  may  be  shown, 
however,  that  a  sufficient  condition  for  a  unique  solution  ...  to  exist 
is  that  a(m)  and  8(m)  not  be  identically  zero  and  that  a(m)  not  be 
related  to  3(M  -  1  -  m)  by  a  constant  scale  factor."  However,  no  proof 
of  that  statement  was  provided.  A  counterexample  to  that  claim  is  shown 
in  Figure  1.  Figures  1(a)  and  1(b)  show  two  different  functions  having 
the  same  boundaries  as  each  other,  and  for  both  objects  a(m)  is  not 
proportional  to  8(M  -  1  -  m),  and  yet  they  have  the  same  Fourier  modulus 
and  the  same  autocorrelation  function,  which  is  shown  in  Figure  1(c). 
Therefore  knowledge  of  the  boundaries  is  not  necessarily  sufficient 
information  for  a  unique  reconstruction. 

An  infinite  number  of  counterexamples  can  be  generated.  From  the 

13 

theory  of  Bruck  and  Sodin  it  is  known  that  the  solution  of  the  phase 

retrieval  problem  [but  not  necessarily  of  Eq.  (5)]  is  unique  unless  the 

Fourier  transform  of  the  object  is  a  factorable  polynomial,  which  is 

unlikely  to  happen  by  chance  for  the  two-dimensional  case. 

Factorability  of  the  Fourier  transform  is  equivalent  to  the  object  being 

expressible  as  a  convolution  of  two  functions,  and  so  ambiguous  cases 

can  be  constructed  by  forming  an  object  by  convolving  (or 

14 

cross-correlating)  two  functions.  The  object  in  Figure  1(a)  was 
fabricated  by  cross-correlating  the  functions  shown  in  Figures  2(a)  and 
2(b).  The  ambiguous  solution  shown  in  Figure  1(b)  is  the  inverted 
convolution  of  the  functions  shown  in  Figures  2(a)  and  2(b).  An 
infinite  number  of  other  ambiguous  examples  can  be  obtained  by  replacing 
the  values  1,  1,  1  and  2  of  the  function  shown  in  Figure  2(b)  by  other 
values,  with  minor  restrictions  on  those  values. 


Terim 


12 

The  recursive  algorithm  involves  the  solution  of  2M  -  1  linear 
equations  in  2M  -  4  unknowns.  One  problem  with  this  is  that  for  m  =  -M 
+  1  and  for  m  =  M  -  1,  Eq.  (5)  involves  only  the  known  boundary  values 
and  not  the  unknowns.  Therefore  one  has  only  2M  -  3  linear  equations  in 
2M  -  4  unknowns  to  begin  with.  A  second  problem  is  that  upon  inspection 
of  those  equations  one  finds  that,  for  the  ambiguous  case  shown  in 
Figure  1,  two  or  more  of  them  are  dependent  equations.  Since  the  number 
of  remaining  linear  independent  equations  is  fewer  than  the  number  of 
unknowns,  the  problem  is  underdetermined  and  multiple  solutions  exist. 

Consider  the  particular  example  of  Figure  1(c),  for  which  one 
searches  for  solutions  of  the  form  shown  in  Figure  2(c),  having  the  a 
priori  known  boundary  values.  Of  the  2M  -  3  =  7  linear  equations  of  Eq. 
(5),  utilizing  the  second  row  of  Figure  1(c),  one  finds  that  three  are 
dependent,  leaving  only  four  independent  equations  in  six  unknowns. 
Therefore  one  can,  for  example,  choose  values  a  and  b  in  Figure  2(c) 
arbitrarily,  and  then  the  values  of  c,  d,  e  and  f  are  determined.  At 
this  point  the  algorithm  of  Reference  12  would  have  been  stopped. 
However,  if  one  continues  to  solve  the  equations  using  the  next  row  of 
the  autocorrelation,  then  one  arrives  at  a  quadratic  equation  in  one  of 
the  variables  yielding  exactly  two  solutions,  those  shown  in  Figures 
1(a)  and  1(b). 

From  the  example  discussed  above  it  is  seen  that  the  recursive 
algorithm  of  Reference  12  is  much  like  the  recursive  algorithm  of 
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Dallas  ,  in  which  a  tree  of  solutions  may  grow  with  each  iteration,  and 
ambiguities  are  resolved  only  if  the  tree  can  be  pruned  in  later 
Iterations. 
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3.  Some  Unique  Cases 

Despite  the  nonuniqueness  demonstrated  in  the  previous  section, 
there  are  some  specific  classes  of  objects  for  which  the  solution  is 
unique.  These  unique  objects  have  supports  (or  shapes)  of  special 
types. 

Certain  classes  of  objects  having  latent  reference  points  can  be 
reconstructed  using  a  simpler  recursive  algorithm  than  the  one  described 
in  the  previous  section.  The  simpler  recursive  algorithm**  selects  the 
order  of  the  equations  being  solved  such  that  at  each  step  one  must 
solve  only  a  single  linear  equation  for  a  single  unknown,  which  is  a 
trivial  computation  that  always  gives  a  unique  result.  It  is  required 
that  no  division  by  zero  be  allowed  and  this  is  ensured  by  the 
requirement  that  the  values  of  the  latent  reference  points  not  be  zero. 
The  latent  reference  points  act  in  a  similar  manner  to  reference  points 
for  holography,  only  they  do  not  initially  satisfy  the  holographic 
separation  condition.  Examples  of  objects  that  can  be  uniquely 
reconstructed  in  this  manner  include  (Fiddy-Brames-Dainty  )  objects 
within  a  rectangle  plus  a  point  off  one  corner  of  the  rectangle,  and 
objects  having  other  supports  as  well.**  In  most  cases  the  support  of 
the  object  must  be  known  a  priori  in  order  to  ensure  that  one  obtains  a 
unique  reconstruction,  since  it  is  usually  not  possible  to  deduce  the 

support  of  the  object  from  the  support  of  its  autocorrelation*^1. 

15 

However,  for  the  Fiddy-Brames-Dainty  objects  the  support  can  be 

deduced  from  the  autocorrelation  support,  and  so  the  reconstruction  in 

that  case  is  unconditionally  unique.**  For  these  cases  the  objects  may 

be  complex-valued,  in  contrast  with  the  restriction  to  real-valued 

objects  for  the  reconstruction  algorithm  discussed  in  the  previous 

section.  Furthermore,  for  these  cases  the  boundary  values  need  not  be 

known  a  priori  since  they  are  computed  in  the  first  step  of  the 

11  12 

recursive  algorithm  ’  . 
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4.  Conclusions 

Although  boundary  conditions  are  a  powerful  constraint  for  the 

phase  retrieval  problem,  it  has  been  proven  by  counterexample  that 

knowledge  of  the  boundary  conditions  (the  values  of  the  edges  of  the 

object)  is  not  sufficient  to  ensure  a  unique  solution.  In  practice  it 

may  be  that  a  unique  solution  is  usually  obtained  simply  because  2-D 

phase  retrieval  is  usually  unique  even  when  the  boundary  conditions  are 
13 

not  known.  It  is  not  yet  known  what  extra  constraints  are  necessary 
to  ensure  uniqueness  in  general. 

What  seems  to  be  more  important  to  ensure  uniqueness  is  that  the 

object's  support  be  a  member  of  a  special  class  of  supports.  It  is  not 

yet  known  in  general  exactly  what  properties  the  support  must  have 

(except  for  the  special  cases  mentioned  in  the  previous  section)  to 

ensure  uniqueness;  but  it  is  known  that  objects  with  separated 

supports*®’*7  are  more  likely  to  be  unique  (even  in  the  1-D  case)  and 

objects  having  complicated  supports  tend  to  be  easier  to  reconstruct 

18 

than  objects  with  convex  symmetric  support  in  the  2-D  case  . 

The  value  of  the  recursive  algorithms  may  be  more  in  their 

predictions  of  uniqueness  than  in  their  ability  to  reconstruct  images, 

11  12 

since  they  tend  to  be  very  sensitive  to  noise  ’.  A  more  stable 
reconstruction  method  would  be  the  iterative  Fourier  transform 
approach®,  which  repeatedly  reinforces  both  the  measured  data  and  the  a 
priori  constraints  on  the  reconstructed  image. 
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FIGURE  CAPTIONS 
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Counterexample  to  uniqueness  claims  .  Two  different  objects 
(a)  and  (b)  have  the  same  boundary  values  and  also  have  the  same 
Fourier  modulus  and  the  same  autocorrelation  (c). 

Functions  (a)  and  (b)  which  generate  the  object  shown  in  Figure 
1(a)  by  cross-correlation  and  in  Figure  1(b)  by  convolution. 

The  general  form  (c)  of  the  objects  which  have  the  autocorrelation 
shown  in  Figure  1(c). 
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